Category Archives: R

八 28

R语言实现并行计算

Posted on by

前面我提到有一个大的运算任务需要很久才完成，所以用到了进度条来监控过程，但并不是改善了计算速度，所以需要用到并行计算，我又在网上找了找。

同样也是一个包，跟matlab的实现过程很像

library(parallel)

cl.cores <- detectCores() #检查当前电脑可用核数。

cl <- makeCluster(cl.cores) #使用刚才检测的核并行运算

#这里用clusterEvalQ或者par开头的apply函数族就可以进行并行计算啦

stopCluster(cl)

R-Doc里这样描述makeCluster函数：Creates a set of copies of R running in parallel and communicating over sockets. 即同时创建数个R进行并行运算。在该函数执行后就已经开始并行运算了，电脑可能会变卡一点。尤其在执行par开头的函数时。

在并行运算环境下，常用的一些计算方法如下：

1、clusterEvalQ(cl,expr)函数利用创建的cl执行expr。这里利用刚才创建的cl核并行运算expr。expr是执行命令的语句，不过如果命令太长的话，一般写到文件里比较好。比如把想执行的命令放在Rcode.r里：clusterEvalQ(cl,source(file="Rcode.r"))

2、par开头的apply函数族。这族函数和apply的用法基本一样，不过要多加一个参数cl。一般如果cl创建如上面cl <- makeCluster(cl.cores)的话，这个参数可以直接用作parApply(cl=cl,…)。当然Apply也可以是Sapply,Lapply等等。注意par后面的第一个字母是要大写的，而一般的apply函数族第一个字母不大写。

另外要注意，即使构建了并行运算的核，不使用parApply()函数，而使用apply()函数的话，则仍然没有实现并行运算。换句话说，makeCluster只是创建了待用的核，而不是并行运算的环境。

参考：http://www.r-bloggers.com/lang/chinese/1131

然后我模仿着用并行计算实现自己的需求

#it did work very fast
library(parallel)
cl.cores <- detectCores()
cl <- makeCluster(cl.cores)
clusterExport(cl, "all_dat_t")  #这里是重点，因为并行计算里面用到了自定义函数
clusterExport(cl, "all_prob_id") #但是这个函数需要用到这两个数据，所以需要把这两个数据加载到并行计算环境里面
prob_202723_s_at=parSapply( #我这里用的parSapply来实现并行计算
cl=cl,  #其中cl是我前面探测到的core数量，

deviation_prob, #deviation_prob是我待并行处理的向量

test_pro #这里其实应该是一个自定义函数，我这里就不写出来了,对上面的deviation_prob向量的每个探针都进行判断
)

八 28

R语言实现进度条

Posted on by

我也是临时在网上搜索到的教程，然后简单看了一下就实现了，其实就是就用到了一个名称为tcltk的包，直接查看函数tkProgressBar就可以知道怎么用啦！

下面是网上的一个小的示例代码（么有实际意义，仅为举例而已）：

library(tcltk2)

u <- 1:2000

plot.new()

pb <- tkProgressBar("进度","已完成 %",  0, 100)

for(i in u) {

x<-rnorm(u)

points(x,x^2,col=i)

info <- sprintf("已完成 %d%%", round(i*100/length(u)))

setTkProgressBar(pb, i*100/length(u), sprintf("进度 (%s)", info), info)

}

close(pb)#关闭进度条

但是下面的代码是我模仿上面这个教程自己实现的。

[R]

# 以下是实现进度条
library(tcltk2)
plot.new()
pb <- tkProgressBar("进度","已完成 %", 0, 100)
prob_202723_s_at_value=rep(0,length(deviation_prob))
start_time=Sys.time() #这里可以计时，因为要实现进度条的一般都是需要很长运算时间
for (i in 1:length(deviation_prob)) {
tmp=test_pro(deviation_prob[i]) #test_pro是我自定义的一个函数，判断该探针是否符合要求。
if (length(tmp)!=0){prob_202723_s_at_value[i]=tmp}
info <- sprintf("已完成 %d%%", round(i*100/length(deviation_prob)))  #进度条就是根据循环里面的i来看看循环到哪一步了
setTkProgressBar(pb, i*100/length(deviation_prob), sprintf("进度 (%s)", info), info)
}
close(pb)#关闭进度条
end_time=Sys.time()
cat(end_time-start_time)

[/R]

八 28

R语言-比较数据框提取列的速度

Posted on by

结论：从数据框里面取某列数据，三种方法的时间消耗区别很大，直接用索引值，是最快的，而用$符号其次，用列名最慢。

我在R里面建立了一个表达量矩阵，列名是一个个样品，行是一个个探针，矩阵值是该探针在该样品测定的表达量。

那么，如果我要看看名为"202723_s_at"的探针的表达向量与其它所有探针的表达向量的相关系数，我可以用以下三种方法：

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t$"202723_s_at",x)))

user  system elapsed

22.93    0.03   23.03

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,"202723_s_at"],x)))

Timing stopped at: 92.02 5.32 97.66

太耗时间了，省去

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,grep(prob,names(all_dat_t))],x)))

Timing stopped at: 13.55 0.04 13.66

> prob_num=grep(prob,names(all_dat_t))

> system.time(apply(all_dat_t,2,function(x)  cor(all_dat_t[,prob_num],x)))

user  system elapsed

8.14    0.01    8.17

可以看出，如果我首先根据探针名，grep出它在该表达量矩阵的列数，然后用列数来提取它的表达量是最快的，而且时间改善非常明显！

我们再探究一下cor函数的效率问题

探究的矩阵有54675个变量，每个变量均有189个观测值，如果取这个大矩阵的部分变量来求相关系数，结果如下！

> system.time(cor(all_dat_t[,1:10]))

user  system elapsed

0.001   0.000   0.001

> system.time(cor(all_dat_t[,1:100]))

user  system elapsed

0.003   0.000   0.003

> system.time(cor(all_dat_t[,1:1000]))

user  system elapsed

0.107   0.002   0.108

> system.time(cor(all_dat_t[,1:10000]))

user  system elapsed

11.102   0.849  11.983

> system.time(cor(all_dat_t)) 约六分钟也是可以搞定的

但是如果cor(all_dat_t)，六分钟后得到的相关系数矩阵约32G，非常恐怖！

但是它很明显没有把这个32G相关系数矩阵存储到内存，因为我的机器本来就16G内存。我至今不能明白R具体实现机理

 

八 27

根据基因表达量对样品进行分类ConsensusClusterPlus

Posted on by 
bioconductor系列的包都是一样的安装方式：
source("http://bioconductor.org/biocLite.R")
biocLite("ConsensusClusterPlus")

这个包是我见过最简单的包， 加载只有做好输入数据，只需要一句话即可运行，然后默认输出所有结果

读这个包的readme，很容易学会
就是做好一个需要来进行分类的样品的表达量矩阵。或者选择上一篇日志用GEOquery这个包下载的表达量矩阵也可以进行分析
因为这个包是用ALL数据来做测试的，所以可以直接加载这个数据结果，这样就能得到表达矩阵啦
library(ALL)
data(ALL)
d=exprs(ALL)
d[1:5,1:5]
可以看到数据集如下

> d[1:5,1:5]

             01005    01010    03002    04006    04007
1000_at   7.597323 7.479445 7.567593 7.384684 7.905312
1001_at   5.046194 4.932537 4.799294 4.922627 4.844565
1002_f_at 3.900466 4.208155 3.886169 4.206798 3.416923
1003_s_at 5.903856 6.169024 5.860459 6.116890 5.687997
1004_at   5.925260 5.912780 5.893209 6.170245 5.615210
> dim(d)
[1] 12625   128
共128个样品，12625个探针数据
也有文献用RNAs-seq的RPKM值矩阵来做
对上面这个芯片表达数据我们一般会简单的进行normalization ，然后取在各个样品差异很大的那些gene或者探针的数据来进行聚类分析
mads=apply(d,1,mad)
d=d[rev(order(mads))[1:5000],]

d = sweep(d,1, apply(d,1,median,na.rm=T))

#也可以对这个d矩阵用DESeq的normalization 进行归一化，取决于具体情况
library(ConsensusClusterPlus)
#title=tempdir() #这里一般改为自己的目录
title="./" #所有的图片以及数据都会输出到这里的
results = ConsensusClusterPlus(d,maxK=6,reps=50,pItem=0.8,pFeature=1,
 title=title,clusterAlg="hc",distance="pearson",seed=1262118388.71279,plot="png")
这样就OK了，你指定的目录下面会输出大于9个图片

clipboard

大家看看说明书就知道这个包的输出文件是什么了。
很多参数都是需要调整的，一般我们的maxK=6是根据实验原理来调整，如果你的样品应该是要分成6类以上，那么你就要把maxK=6调到一点。
查看结果results[[2]][["consensusClass"] 可以看到各个样品被分到了哪个类别里面去
results[[3]][["consensusClass"]
results[[4]][["consensusClass"] 等等
八 27

从GEO数据库下载矩阵数据-可以直接进行下游分析

Posted on by 
bioconductor系列的包都是一样的安装方式：
source("http://bioconductor.org/biocLite.R")
    biocLite("GEOquery")

以前GEO数据库主要是microarray的芯片数据，后来有了RNA-seq，如果同时做多个样品的RNA-seq，表达量矩阵后来也可以上传到GEO数据库里面，只有看到文献里面有提到GEO数据库，都可以通过这个R包俩进行批量下载，其实就是网页版的一个API调用而已：

GEO数据库里面有四种数据

At the most basic level of organization of GEO, there are four basic entity types.
 The first three (Sample, Platform, and Series) are supplied by users;
the fourth, the dataset, is compiled and curated by GEO sta from the user-submitted data.
GEO accession number (GPLxxx). 
GEO accession number (GSMxxx)
GEO accession number (GSExxx). 
GEO DataSets (GDSxxx)
记住大小关系：一个GDS可以有多个GSM，一个GSM可以有多个GSE，至于GPL，一般不接触的
我们通常接触的都是GSE系列（一个GSE里面有多个GSM）的数据，而且这个包最重要的就是一个getGEO函数。
只要你通过文献确定了你的检索号，就可以通过这个函数来下载啦
检索号一般是A character string representing a GEO object for download and
          parsing.  (eg., 'GDS505','GSE2','GSM2','GPL96'

这个函数有很多参数，除非你需要下载的文件，那么就设置destdir到你喜欢的目录，如果只需要表达量数据就不用了。

 getGEO(GEO = NULL, filename = NULL, destdir = tempdir(), GSElimits=NULL,
     GSEMatrix=TRUE,AnnotGPL=FALSE)
例如：
gds <- getGEO("GDS10") 会返回一个对象，而且下载数据一般会在tmp目录下面，当然如果你需要保存这些文件，你可以自己制定下载目录及文件名。
gse2553 <- getGEO("GSE2553")
GDS2eSet函数可以把上面这个下载函数得到的对象(要确定是GDS而不是GSE)变成表达对象
pData和exprs函数都可以处理上面这个表达对象，从而分别得到样品描述矩阵和样品表达量矩阵
综合一起就是
g4100 <- GDS2eSet(getGEO("GDS4100"))
g4102 <- GDS2eSet(getGEO("GDS4102"))
e4102<-exprs(g4102)
e4100<-exprs(g4100)
这样的代码，这个e4100和e4102就都是一个数值矩阵啦，可以进行下游分析，但是如果是下载的GSM数据
就用下面这个代码，GSE26253_series_matrix.txt是通过GSEMatrix=TRUE这个参数特意下载到你的目录的
expr_dat=read.table("GSE26253_series_matrix.txt",comment.char="!",stringsAsFactors=F)
这样读取也是一个数值矩阵
具体大家可以看这个包的说明书
#Download GDS file, put it in the current directory, and load it:
gds858 <- getGEO('GDS858', destdir=".")
如果使用了GSEMatrix=TRUE这个参数，那么除了下载soft文件，还有表达量矩阵文件，可以直接用read.table读取那个文件。
#Or, open an existing GDS file (even if its compressed):
gds858 <- getGEO(filename='GDS858.soft.gz')
下面这个下载的是GSE对象，GDS对象还有大一点
GEOquery-下载结果gds对象
七 21

R语言用hclust进行聚类分析

Posted on by

聚类的基础就是算出所有元素两两间的距离，我们首先做一些示例数据，如下：

x=runif(10)

y=runif(10)

S=cbind(x,y)                                 #得到2维的数组

rownames(S)=paste("Name",1:10,"")             #赋予名称，便于识别分类

out.dist=dist(S,method="euclidean")           #数值变距离

这个代码运行得到的S是一个矩阵，如下

> S

x         y

Name 1   0.41517985 0.4697017

Name 2   0.35653781 0.1132367

Name 3   0.52253349 0.3680286

Name 4   0.80558684 0.9834687

Name 5   0.04564145 0.8560690

Name 6   0.11044397 0.2988598

Name 7   0.34984447 0.8515141

Name 8   0.28097709 0.1260050

Name 9   0.81771888 0.5976135

Name 10 0.40700158 0.5236567

可以看出里面共有10个点，它们的X,Y坐标均已知，我们有6总方法可以求矩阵

image001

注释：在聚类中求两点的距离有：

1，绝对距离：manhattan

2，欧氏距离：euclidean 默认

3，闵科夫斯基距离：minkowski

4，切比雪夫距离：chebyshev

5，马氏距离：mahalanobis

6，蓝氏距离：canberra

用默认的算法求出距离如下

算出距离后就可以进行聚类啦！

out.hclust=hclust(out.dist,method="complete") #根据距离聚类

注释：聚类也有多种方法：

1，类平均法：average

2，重心法：centroid

3，中间距离法:median

4，最长距离法：complete 默认

5，最短距离法：single

6，离差平方和法：ward

7，密度估计法：density

接下来把聚类的结果图画出来

plclust(out.hclust)                           #对结果画图

image003

rect.hclust(out.hclust,k=3)                   #用矩形画出分为3类的区域

image005

out.id=cutree(out.hclust,k=3)                 #得到分为3类的数值

这里的out.id就是把每个点都分类了的分类数组，1,2,3.

 

七 21

用R语言批量做T检验

Posted on by

需要做T检验的的数据如下：其中加粗加黑的是case，其余的是control，需要进行检验的是

p_ma    p_mg    p_ag    p_mag 这四列数据，在case和control里面的差异，当然用肉眼很容易看出，control要比case高很多

 

individual m a g ma mg ag mag p_ma p_mg p_ag p_mag
HBV10-1 33146 3606 2208 308 111 97 39 0.79% 0.29% 0.25% 0.10%
HBV15-1 23580 10300 3140 1469 598 560 323 4.19% 1.71% 1.60% 0.92%
HBV3-1 107856 26445 15077 4773 2383 1869 1130 3.34% 1.67% 1.31% 0.79%
HBV4-1 32763 6448 4384 579 291 295 124 1.35% 0.68% 0.69% 0.29%
HBV6-1 122396 6108 7953 911 796 347 144 0.67% 0.59% 0.26% 0.11%
IGA-1 31337 22167 14195 3851 2752 4101 2028 6.50% 4.65% 6.92% 3.42%
IGA-10 6713 9088 12801 2275 2470 4284 1977 10.54% 11.44% 19.85% 9.16%
IGA-11 61574 23622 15076 5978 4319 3908 2618 6.71% 4.84% 4.38% 2.94%
IGA-12 38510 23353 20148 6941 6201 6510 4328 10.38% 9.27% 9.73% 6.47%
IgA-13 155379 81980 65315 26055 20085 17070 10043 10.38% 8.00% 6.80% 4.00%
IgA-14 345430 86462 71099 27541 21254 10923 6435 6.06% 4.67% 2.40% 1.42%
IgA-15 3864 3076 1942 378 207 389 145 4.66% 2.55% 4.80% 1.79%
IgA-16 3591 4106 2358 427 174 424 114 4.64% 1.89% 4.61% 1.24%
IgA-17 893 1442 799 68 28 78 18 2.27% 0.94% 2.61% 0.60%
IGA-2 23097 5083 5689 910 951 1173 549 2.89% 3.02% 3.72% 1.74%
IGA-3 14058 9364 8565 2130 1953 2931 1436 8.03% 7.36% 11.05% 5.41%
IGA-4 81571 36894 33664 11346 10131 9908 6851 8.86% 7.91% 7.74% 5.35%
IGA-5 27626 6492 4503 963 752 942 410 2.64% 2.06% 2.58% 1.12%
IGA-7 55348 5956 4028 833 476 468 207 1.30% 0.74% 0.73% 0.32%
IGA-8 31671 17097 10443 3894 2666 3514 2003 7.56% 5.17% 6.82% 3.89%

 

但是我们需要写程序对这些百分比在case和control里面进行T检验。

a=read.table("venn_dat.txt",header=T)

type=c(rep("control",5),rep("case",12),"control","control","case")

for (i in 9:12)

{

dat=as.numeric(unlist(lapply(a[i],function(x) strsplit(as.character(x),"%"))))

filename=paste(names(a)[i],".png",sep="")

png(file=filename, width = 320, height = 360)

b=t.test(dat~type)

txt=paste("p-value=",round(b$p.value[1],digits=4),sep="")

plot(as.factor(type),dat,ylab="percent(%)",main=names(a)[i],sub=txt,cex.axis=1.5,cex.sub=2,cex.main=2,cex.lab=1.5)

dev.off()

}

得到的图片如下：

image008

七 21

用R语言批量画韦恩图

Posted on by

需要画韦恩图的文件如下所示：

#CDR3_aa    count_all    count_IgM    count_IgA    count_IgG
CARGVDAGVDYW    243    25    196    22
CARHPRNYGNFDYW    174    171    3    0
CARENTMVRGVINPLDYW    166    8    75    83
CAREASDSISNWDDWYFDLW    129    15    114    0
CARDPDNSGAFDPW    118    1    117    0
CAKDLGGYW    98    3    4    91
CAREVADYDTYGWFLDLW    95    26    68    1
CVRNRGFFGLDIW    82    0    1    81
CARRSTNYHGWDYW    80    3    2    74

此处省略一万行。

简单解释一下数据，第一列是CDR3序列，我们需要对count_IgM    count_IgA    count_IgG这三列数据进行画韦恩图，数字大于0代表有，数字为0代表无。

这样我们根据序列就能得出每列数据所有的CDR3序列，即不为0的CDR3序列

每个个体都会输出一个统计文件，共20个文件需要画韦恩图

image005

对这个统计文件就可以进行画韦恩图。

画韦恩图的R代码如下：

library(VennDiagram)

files=list.files(path = ".", pattern = "type")

for (i in files){

a=read.table(i)

individual=strsplit(i,"\\.")[[1]][1]

image_name=paste(individual,".tiff",sep="")

IGM=which(a[,3]>0)

IGA=which(a[,4]>0)

IGG=which(a[,5]>0)

venn.diagram(list(IGM=IGM,IGA=IGA,IGG=IGG), fill=c("red","green","blue"), alpha=c(0.5,0.5,0.5), cex=2, cat.fontface=4, fontfamily=3, filename=image_name)

}

image007

但事实上，这个韦恩图很难表现出什么，因为我们的每个个体的count_IgM    count_IgA    count_IgG总数不一样。

七 19

Bioconductor系列之GenomicAlignments

Posted on by

Bioconductor系列包的安装方法都一样

source("http://bioconductor.org/biocLite.R")biocLite(“GenomicAlignments”)

这个包设计就是用来处理bam格式的比对结果的，具体功能非常多，其实还不如自己写脚本来做这些工作，更方便一点，实在没必要学别人的语法。大家有兴趣可以看GenomicAlignments的主页介绍，三个pdf文档来介绍。

http://www.bioconductor.org/packages/release/bioc/html/GenomicAlignments.html

我们首先读取示例pasillaBamSubset包带有的bam文件

library(pasillaBamSubset)

un1 <- untreated1_chr4()  # single-end reads

library(GenomicAlignments)

reads1 <- readGAlignments(un1)

查看reads1这个结果，可以看到把这个bam文件都读成了一个数据对象GAlignments object，

Bioconductor系列之GenomicAlignments632

readGAlignments这个函数就是读取我们bam文件的，生成的对象可以用多个函数来继续查看信息，比如它的列名都是函数

Seqnames(),strand(),cigar(),qwidth(),start(),end(),width(),njunc() 这些函数对这个GAlignments对象处理后会得到比对的各个情况的向量。

我们也可以读取这个GenomicAlignments包自带的bam文件，而不是pasillaBamSubset包带有的bam文件

> fls <- list.files(system.file("extdata", package="GenomicAlignments"),+                   recursive=TRUE, pattern="*bam$", full=TRUE)

> fls

[1] "C:/Program Files/R/R-3.1.2/library/GenomicAlignments/extdata/sm_treated1.bam"  [2] "C:/Program Files/R/R-3.1.2/library/GenomicAlignments/extdata/sm_untreated1.bam"

reads1 <- readGAlignments(fls[1])

比如我随意截图 cigar(reads1)的结果，就可以看出，大部分都是75M这样的比对情况，可以对这个向量进行统计完全匹配的情况，总共有6077总比对情况。

unique(cigar(reads1))

Bioconductor系列之GenomicAlignments1315

也可以用table格式来查看cigar，可以看到大部分都是M

head(cigarOpTable(cigar(reads1)))

Bioconductor系列之GenomicAlignments1385

接下来我们再读取PE测序数据的比对bam结果，看看分析方法有哪些。

Bioconductor系列之GenomicAlignments1422

U3.galp <- readGAlignmentPairs(untreated3_chr4(), use.names=TRUE, param=param0)head(U3.galp)

也可以分开查看左右两端测序数据的比对情况，跟上面的head是对应关系，上面的

[169,  205] -- [ 326,  362]在下面就被分成左右两个比对情况了。

> head(first(U3.galp))

GAlignments object with 6 alignments and 0 metadata columns:      seqnames strand       cigar    qwidth     start       end     width     njunc         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>  [1]     chr4      +         37M        37       169       205        37         0  [2]     chr4      +         37M        37       943       979        37         0  [3]     chr4      +         37M        37       944       980        37         0  [4]     chr4      +         37M        37       946       982        37         0  [5]     chr4      +         37M        37       966      1002        37         0  [6]     chr4      +         37M        37       966      1002        37         0  -------  seqinfo: 8 sequences from an unspecified genome

> head(last(U3.galp))

GAlignments object with 6 alignments and 0 metadata columns:      seqnames strand       cigar    qwidth     start       end     width     njunc         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer> <integer>  [1]     chr4      -         37M        37       326       362        37         0  [2]     chr4      -         37M        37      1086      1122        37         0  [3]     chr4      -         37M        37      1119      1155        37         0  [4]     chr4      -         37M        37       986      1022        37         0  [5]     chr4      -         37M        37      1108      1144        37         0  [6]     chr4      -         37M        37      1114      1150        37         0  -------  seqinfo: 8 sequences from an unspecified genome

用isProperPair可以查看pe测序数据比对情况的，哪些没有配对比对成功

> table(isProperPair(U3.galp))FALSE  TRUE 29518 45828

还可以用Rsamtools包对我们的bam文件进行统计，看下面的代码中fls[1]和un1都是上面得到的bam文件全路径。

> library(Rsamtools)

> quickBamFlagSummary(fls[1], main.groups.only=TRUE)

group |    nb of |    nb of | mean / max                                   of |  records |   unique | records per                              records | in group |   QNAMEs | unique QNAMEAll records........................ A |     1800 |     1800 |    1 / 1  o template has single segment.... S |     1800 |     1800 |    1 / 1  o template has multiple segments. M |        0 |        0 |   NA / NA      - first segment.............. F |        0 |        0 |   NA / NA      - last segment............... L |        0 |        0 |   NA / NA      - other segment.............. O |        0 |        0 |   NA / NA Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.Indentation reflects this.

> library(Rsamtools)

> quickBamFlagSummary(un1, main.groups.only=TRUE)

group |    nb of |    nb of | mean / max                                   of |  records |   unique | records per                              records | in group |   QNAMEs | unique QNAMEAll records........................ A |   204355 |   190770 | 1.07 / 10  o template has single segment.... S |   204355 |   190770 | 1.07 / 10  o template has multiple segments. M |        0 |        0 |   NA / NA      - first segment.............. F |        0 |        0 |   NA / NA      - last segment............... L |        0 |        0 |   NA / NA      - other segment.............. O |        0 |        0 |   NA / NA Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.Indentation reflects this.

 

下面讲最后一个综合应用

biocLite("RNAseqData.HNRNPC.bam.chr14")

#这是一个665.5MB的RNA-seq测序数据的比对结果

library(RNAseqData.HNRNPC.bam.chr14)bamfiles <- RNAseqData.HNRNPC.bam.chr14_BAMFILESnames(bamfiles)  # the names of the runslibrary(Rsamtools)quickBamFlagSummary(bamfiles[1], main.groups.only=TRUE)flag1 <- scanBamFlag(isFirstMateRead=TRUE, isSecondMateRead=FALSE,                     isDuplicate=FALSE, isNotPassingQualityControls=FALSE)param1 <- ScanBamParam(flag=flag1, what="seq") #这里是设置readGAlignments的读取参数library(GenomicAlignments)gal1 <- readGAlignments(bamfiles[1], use.names=TRUE, param=param1)

mcols(gal1)$seqoqseq1 <- mcols(gal1)$seqis_on_minus <- as.logical(strand(gal1) == "-")oqseq1[is_on_minus] <- reverseComplement(oqseq1[is_on_minus])

这样就还原得到了原始测序数据啦！

七 19

Bioconductor系列之pasillaBamSubset

Posted on by

这个包主要有bam文件测试数据
> biocLite("pasillaBamSubset")
BioC_mirror: http://bioconductor.orgUsing Bioconductor version 3.0 (BiocInstaller 1.16.5), R version 3.1.2.
Installing package(s) 'pasillaBamSubset'
trying URL 'http://bioconductor.org/packages/3.0/data/experiment/bin/windows/contrib/3.1/pasillaBamSubset_0.3.1.zip'
Content type 'application/zip' length 31514402 bytes (30.1 Mb)
打开pasillaBamSubset包的安装地址就可以看到里面有几个bam文件
Several functions are available for reading BAM files into R:
而且加载包的同时也引入了几个读取bam文件的函数
readGAlignments()
readGAlignmentPairs()
readGAlignmentsList()
scanBam()
Bioconductor系列之pasillaBamSubset448
加载包就可以看到用两个函数得到包自带的数据文件的地址，主要是有很多人不一定把包安装在C盘，所以用函数来定位文件更加安全一点
> library(pasillaBamSubset)
> untreated1_chr4()
[1] "C:/Program Files/R/R-3.1.2/library/pasillaBamSubset/extdata/untreated1_chr4.bam"
> untreated3_chr4()
[1] "C:/Program Files/R/R-3.1.2/library/pasillaBamSubset/extdata/untreated3_chr4.bam"

接下来我们就看看如何读取这些bam文件的
library(pasillaBamSubset)
un1 <- untreated1_chr4() # single-end reads library(GenomicAlignments) reads1 <- readGAlignments(un1) cvg1 <- coverage(reads1) 查看reads1这个结果，可以看到把这个bam文件都读成了一个数据对象GAlignments object， Bioconductor系列之pasillaBamSubset1142
针对着个数据对象有很多操作，其中一个coverage操作是来自于GenomicFeatures
或者GenomicAlignments函数的，可以算出测序覆盖情况。
可以看到这个bam文件里面的比对情况大多几种在4号染色体里面
> cvg1$chr4
integer-Rle of length 1351857 with 122061 runs
Lengths: 891 27 5 12 13 45 5 12 13 ... 5 1 1 3 10
Values : 0 1 2 3 4 5 4 3 2 ... 12 11 10 6
> mean(cvg1$chr4)
[1] 11.33746
> max(cvg1$chr4)[1] 5627
可以看到平均测序深度是11.3X，最大测序深度是5627X

七 18

Bioconductor系列之biomaRt

Posted on by

biomaRt是一个超级网络资源库，里面的信息非常之多，就是网页版的biomaRt的R语言接口。谷歌搜索 the biomart user’s guide filetype:pdf 这个关键词，就看到关于这个包的详细介绍以及例子，我这里简单总结一下它的用法。

这个包一直在更新，函数总是变化，而且需要联网才能使用，基本上等于废物一个，希望大家不要使用它~

包的安装

Continue reading →

七 18

Bioconductor系列之hgu95av2.db

Posted on by

这个包主要都是芯片数据处理方面的应用，本来我是懒得看的，但是我无意中浏览了一下readme，发现它挺适合被作为数据库的典范来学习。
source("http://bioconductor.org/biocLite.R")
biocLite("org.Hs.eg.db")
biocLite("hgu95av2.db")
还是老办法安装了hgu95av2.db之后用ls可以查看这个包里面有36个映射数据，都是把芯片探针ID号转为其它36种主流ID号的映射而已。
library(hgu95av2.db)
> ls("package:hgu95av2.db")
[1] "hgu95av2" "hgu95av2ACCNUM" "hgu95av2ALIAS2PROBE"
[4] "hgu95av2CHR" "hgu95av2CHRLENGTHS" "hgu95av2CHRLOC"
[7] "hgu95av2CHRLOCEND" "hgu95av2.db" "hgu95av2_dbconn"
[10] "hgu95av2_dbfile" "hgu95av2_dbInfo" "hgu95av2_dbschema"
[13] "hgu95av2ENSEMBL" "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"
[16] "hgu95av2ENZYME" "hgu95av2ENZYME2PROBE" "hgu95av2GENENAME"
[19] "hgu95av2GO" "hgu95av2GO2ALLPROBES" "hgu95av2GO2PROBE"
[22] "hgu95av2MAP" "hgu95av2MAPCOUNTS" "hgu95av2OMIM"
[25] "hgu95av2ORGANISM" "hgu95av2ORGPKG" "hgu95av2PATH"
[28] "hgu95av2PATH2PROBE" "hgu95av2PFAM" "hgu95av2PMID"
[31] "hgu95av2PMID2PROBE" "hgu95av2PROSITE" "hgu95av2REFSEQ"
[34] "hgu95av2SYMBOL" "hgu95av2UNIGENE" "hgu95av2UNIPROT"
但是用这个包自带的函数capture.output(hgu95av2())可以看到这些映射并没有囊括我们标准的hg19版本的2.3万个基因，也就说这个芯片设计的探针只有1.1万个左右。但它根据已有的org.Hs.eg.db来构建了它自己芯片独有的数据包，这样就显得更加正式规范化了，这启发我们研发人员在做一些市场产品的同时也可以把自己的研究成果通主流数据库数据形式结合起来，更易让他人接受。
> capture.output(hgu95av2())
[1] "Quality control information for hgu95av2:"
[2] ""
[3] ""
[4] "This package has the following mappings:"
[5] ""
[6] "hgu95av2ACCNUM has 12625 mapped keys (of 12625 keys)"
[7] "hgu95av2ALIAS2PROBE has 33755 mapped keys (of 103735 keys)"
[8] "hgu95av2CHR has 11540 mapped keys (of 12625 keys)"
[9] "hgu95av2CHRLENGTHS has 93 mapped keys (of 93 keys)"
[10] "hgu95av2CHRLOC has 11474 mapped keys (of 12625 keys)"
[11] "hgu95av2CHRLOCEND has 11474 mapped keys (of 12625 keys)"
[12] "hgu95av2ENSEMBL has 11460 mapped keys (of 12625 keys)"
[13] "hgu95av2ENSEMBL2PROBE has 9814 mapped keys (of 28553 keys)"
[14] "hgu95av2ENTREZID has 11543 mapped keys (of 12625 keys)"
[15] "hgu95av2ENZYME has 2125 mapped keys (of 12625 keys)"
[16] "hgu95av2ENZYME2PROBE has 786 mapped keys (of 975 keys)"
[17] "hgu95av2GENENAME has 11543 mapped keys (of 12625 keys)"
[18] "hgu95av2GO has 11245 mapped keys (of 12625 keys)"
[19] "hgu95av2GO2ALLPROBES has 17214 mapped keys (of 18826 keys)"
[20] "hgu95av2GO2PROBE has 12920 mapped keys (of 14714 keys)"
[21] "hgu95av2MAP has 11519 mapped keys (of 12625 keys)"
[22] "hgu95av2OMIM has 10541 mapped keys (of 12625 keys)"
[23] "hgu95av2PATH has 5374 mapped keys (of 12625 keys)"
[24] "hgu95av2PATH2PROBE has 228 mapped keys (of 229 keys)"
[25] "hgu95av2PMID has 11529 mapped keys (of 12625 keys)"
[26] "hgu95av2PMID2PROBE has 372382 mapped keys (of 432400 keys)"
[27] "hgu95av2REFSEQ has 11506 mapped keys (of 12625 keys)"
[28] "hgu95av2SYMBOL has 11543 mapped keys (of 12625 keys)"
[29] "hgu95av2UNIGENE has 11533 mapped keys (of 12625 keys)"
[30] "hgu95av2UNIPROT has 11315 mapped keys (of 12625 keys)"
[31] ""
[32] ""
[33] "Additional Information about this package:"
[34] ""
[35] "DB schema: HUMANCHIP_DB"
[36] "DB schema version: 2.1"
[37] "Organism: Homo sapiens"
[38] "Date for NCBI data: 2014-Sep19"
[39] "Date for GO data: 20140913"
[40] "Date for KEGG data: 2011-Mar15"
[41] "Date for Golden Path data: 2010-Mar22"
[42] "Date for Ensembl data: 2014-Aug6"
这个hgu95av2.db所加载的36个包都是一种特殊的对象，但是可以把它当做list来操作，是一种类似于hash的对应表格，其中keys是独特的，但是value可以有多个。
既然是类似于list，那我就简单讲几个小技巧来操作这些数据对象。所有的操作都要用as.list()函数来把数据展现出来
> as.list(hgu95av2ENZYME[1])
$`1000_at`
[1] "2.7.11.24"
可以看到这样就提取出来了hgu95av2ENZYME的第一个元素，key是`1000_at`，它所映射的value是 "2.7.11.24"这个酶。
同理对list取元素的三个操作在这里都可以用
> as.list(hgu95av2ENZYME['1000_at'])
$`1000_at`
[1] "2.7.11.24"
> as.list(hgu95av2ENZYME$'1000_at')
[[1]]
[1] "2.7.11.24"
然而这只是list的操作，我们还有一个函数专门对这个对象来取元素，那就是get函数，取多个元素用mget，类似于as.list(hgu95av2ENZYME[1:10])
用函数mget(probes_id,hgu95av2ENZYME)就可以根据自己的probes_id这个向量来选取数据了，当然也要用as.list()来展示出来。
值得一提的是这里有个非常重要的函数是revmap，可以把我们的key和value调换位置。
因为我们的数据对象里面的对应关系不是一对一，而是一对多，例如一个基因往往有多个go通路信息，例如这个就有十几个
as.list(hgu95av2GO$'1000_at')
$`GO:0000165`
$`GO:0000165`$GOID
[1] "GO:0000165"

$`GO:0000165`$Evidence
[1] "TAS"

$`GO:0000165`$Ontology
[1] "BP"

$`GO:0000186`
$`GO:0000186`$GOID
[1] "GO:0000186"

$`GO:0000186`$Evidence
[1] "TAS"

$`GO:0000186`$Ontology
[1] "BP"
一层层的数据结构非常清晰，但是数据太多，所以它给了一个toTable函数来格式化，就是把一对多的list压缩成表格。
> toTable(hgu95av2GO[1])
probe_id go_id Evidence Ontology
1 1000_at GO:0000165 TAS BP
2 1000_at GO:0000186 TAS BP
3 1000_at GO:0000187 TAS BP
4 1000_at GO:0002224 TAS BP
5 1000_at GO:0002755 TAS BP
6 1000_at GO:0002756 TAS BP
7 1000_at GO:0006360 TAS BP
------------------------------------------------------------
81 1000_at GO:0004707 NAS MF
82 1000_at GO:0001784 IEA MF
83 1000_at GO:0005524 IEA MF
这样就非常方便的查看啦。
当然对这个数据对象还有很多实用的操作。length(),Llength(),Rlength(),keys(),Lkeys(),Rkeys()等等，还有count.mappedkeys(),mappedLkeys()，还有个比较特殊的isNA来判断是否有些keys探针没有对应任何信息。
而且以上这些函数不仅可以用来获取数据对象的信息，还可以修改这个数据对象。
Lkeys(hgu95av2CHR) 可以查看这个对象有12625个探针。
而> Rkeys(hgu95av2CHR)
[1] "19" "12" "8" "14" "3" "2" "17" "16" "9" "X" "6" "1" "7" "10" "11"
[16] "22" "5" "18" "15" "Y" "20" "21" "4" "13" "MT" "Un"
可以看到这些探针分布在不同的染色体上面，如果我这时候给
x=hgu95av2CHR
Rkeys(x) = c("1","2")
toTable(x)可以看到这时候x只剩下1946个探针了，就是1,2号染色体上面的基因了。
然后可以一个个看这些函数的用法，其实就是SQL的增删改查的基本操作而已。
> length(x)
[1] 12625
> Llength(x)
[1] 12625
> Rlength(x)
[1] 2
> head(keys(x))
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at" "1005_at"
> head(Lkeys(x))
[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at" "1005_at"
> head(Rkeys(x))
[1] "1" "2"
> count
count.fields count.mappedLkeys countOverlaps counts<- count.links count.mappedRkeys countQueryHits countSubjectHits count.mappedkeys countMatches counts > count.mappedkeys(x)
[1] 1946
> count.mappedLkeys(x)
[1] 1946
> count.mappedRkeys(x)
[1] 2

七 17

转-R语言内存管理

Posted on by

R中的对象（比如矩阵）在内存中存于两种不同的地方:

第一种是堆内存（heap），其基本单元是“Vcells”，每个大小为8字节，新来一个对象就会申请一块空间，把值全部存在这里，和C里面的堆内存很像;

第二种是地址对（cons cells），主要用来存储地址信息，最小单元一般在32位系统中是28字节、64位系统中是56字节。

1、ls()来查看当前所有对象名，对于每一个对象，可以通过object.size(x)来查看其占用内存的大小。

如果是因为当前对象占用内存过多，那么可以通过处理对象来获取更大的可用内存。一个很有用的方法是改变对象的存储模式，通过storage.mode（x）可以看到某个对象的存储模式，比如某个矩阵默认就是“double”的，如果这个矩阵的数值都是整数甚至0-1，完全没必要使用double来占用空间，可以使用storage.mode(x） <- "integer"将其改为整数型，可以看到该对象的大小会变为原来的一半。

2、object.size()看每个变量占多大内存。

3、memory.size()查看现在的work space的内存使用

4、memory.limit()查看系统规定的内存使用上限。如果现在的内存上限不够用，可以通过memory.limit(newLimit)更改到一个新的上限。注意，在32位的R中，封顶上限为4G，无法在一个程序上使用超过4G （数位上限）。这种时候，可以考虑使用64位的版本。

 

对于一些很大的但无用的中间变量，养成清理的习惯：

可以使用rm(object)删除变量，但是记住，rm后记得使用gc()做Garbage collection，否则内存是不会自动释放的，相当于你没做rm.

七 16

模拟测序lambda_virus基因组

Posted on by
lambda_virus基因组文件是bowtie软件自带的测试数据，共48502个bp，首先我用脚本模拟出它的全打断文件！

perl -alne '{next if /^>/;$a.=$_;}END{$len=length $a;print substr($a.$a,$_,120) foreach 0..$len}' lambda_virus.fa >lamb_virus.120bp

长度均为120bp的片段。

我测序的策略是CTAG碱基重复30次，共加入120个碱基。

对每个120bp片段来说，如果遇到互补碱基就加上，直到120个碱基加完，这样如果比较巧合的话，会有部分碱基能全部加满120bp的，但是如果每个120bp片段的ATCG分布均匀，那么就都应该30bp碱基能被加上。

image001

[perl]
while (<>) {

$seq=$_;$sum=0;

foreach $i (0..120){

$str=substr($seq,$i,2);

if ($str eq "GG"| $str eq "CC"| $str eq "AA"| $str eq "TT"){$sum+=4;}

elsif ($str eq "GT"| $str eq "CG"| $str eq "AC"| $str eq "TA"){$sum+=3;}

elsif ($str eq "GA"|$str eq "CT"| $str eq "AG"| $str eq "TC"){$sum+=2;}

else{$sum+=1;};

#print "$sum\n";

if ($sum>120){print "$i\n";last;}

}

}

[/perl]

perl length.pl lambda_virus.120bp >length.txt

得到结果如下：

 

Length 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
Count 2 19 34 110 204 432 878 1495 2237 3202 4343 5179 5697 5429 4865 4214
Length 52 53 54 55 56 57 58 59 60 61 62 63 64
Count 3249 2499 1735 1090 657 396 228 141 90 48 18 9 3

右表可以看出，大部分测序得到碱基长度都集中在46bp到51bp之间

画出箱线图如下

image003

画出条形图如下：

image005

 

然后我模拟了一个6000bp的基因组，做同样的模拟测序看看评价测序长度分布情况：

Length 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Count 9 22 96 207 322 382 479 671 770 714 706 546 424 232 182 100 52 30 14 21
Length 59 60 61                                  
Count 15 5 2                                  

可以看出这次的测序片段集中在45到51bp

五 29

用R语言的RCurl包结合XML包批量下载生信课件

Posted on by

首先是宾夕法尼亚州立大学（The Pennsylvania State University缩写PSU)的生信课件下载，这个生信不仅有课件，而且在中国的优酷视频网站里面还有全套授课视频，非常棒！

image001

课程主页是http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html

可以看出所有的课件pdf链接都在这一个页面，所以是非常简单的代码！

下面是R代码：

library(XML)

library(RCurl)

library(dplyr)

psu_edu_url='http://www.personal.psu.edu/iua1/courses/2013-BMMB-597D.html';

wp=getURL(psu_edu_url)

base='http://www.personal.psu.edu/iua1/courses/file';

#pse_edu_links=getHTMLLinks(psu_edu_url)

psu_edu_links=getHTMLLinks(wp)

psu_edu_pdf=psu_edu_links[grepl(".pdf$",psu_edu_links,perl=T)]

for (pdf in psu_edu_pdf){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

因为这三十个课件都是接近于10M，所以下载还是蛮耗时间的

image003

其实R语言里面有这个down_url函数，可以直接下载download.file(down_url,filename)

然后我开始下载德国自由大学的生信课件，这次不同于宾夕法尼亚州立大学的区别是，课程主页里面是各个课题的链接，而pdf讲义在各个课题里面，所以我把pdf下载写成了一个函数对我们的课题进行批量处理

library(XML)

library(RCurl)

library(dplyr)

base="http://www.mi.fu-berlin.de/w/ABI/Genomics12";

down_pdf=function(url){

links=getHTMLLinks(url)

pdf_links=links[grepl(".pdf$",links,perl=T)]

for (pdf in pdf_links){

down_url=getRelativeURL(pdf,base)

filename=last(strsplit(pdf,"/")[[1]])

cat("Now we down the ",filename,"\n")

#pdf_file=getBinaryURL(down_url)

#FH=file(filename,"wb")

#writeBin(pdf_file,FH)

#close(FH)

download.file(down_url,filename)

}

}

down_pdf(base)

list_lecture= paste("http://www.mi.fu-berlin.de/w/ABI/GenomicsLecture",1:15,"Materials",sep="")

for ( url in list_lecture ){

cat("Now we process the ",url ,"\n")

try(down_pdf(url))

}

image005

同样也是很多pdf需要下载

接下来下载Minnesota大学的关于生物信息的教程的ppt合集

主页是： https://www.msi.umn.edu/tutorial-materials

 

这个网页里面有64篇pdf格式的ppt，还有几个压缩包，本来是准备写爬虫来爬去的，但是后来想了想有点麻烦，而且还不一定会看，反正也是玩玩

就用linux的命令行简单实现了这个爬虫功能。

curl https://www.msi.umn.edu/tutorial-materials >tmp.txt

perl -alne '{/(https.*?pdf)/;print $1 if $1}' tmp.txt >pdf.address

perl -alne '{/(https.*?txt)/;print $1 if $1}' tmp.txt

perl -alne '{/(https.*?zip)/;print $1 if $1}' tmp.txt >zip.address

wget -i pdf.address

wget -i pdf.zip

这样就可以啦！

 

用爬虫也就是几句话的事情，因为我已经写好啦下载函数，只需要换一个主页即可下载页面所有的pdf文件啦！

 

五 29

R语言对Ozone数据处理笔记

Posted on by

一．首先加载一些包，这样才能获得该ozone数据

数据集介绍：

Ozone数据集是一个三维数组，记录了24×24个空间网格内，从1995年1月到2000年12月，共72个时间点上，中美洲每月的平均臭氧水平。

前两维分别表示纬度和经度，第三维表示时间。

加载包的代码如下：

library("MASS", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("ggplot2", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("plyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("dplyr", lib.loc="C:/Program Files/R/R-3.1.1/library")

library("reshape2", lib.loc="C:/Program Files/R/R-3.1.1/library")

二．我们首先简单看看第一个地点的72个月的臭氧水平变化图。

plot(1:72,ozone[1,1,],type="l")

box(lty = '1373', col = 'red')

grid(nx=NA,ny=NULL,lty=1,lwd=1,col="gray")

看起来还算是蛮有规律的。

image001

三．然后把这72个月的数据分成年份来画图

绘图第一种方式如下：

value <-ozone[1, 1, ]

plot(value[1:12],type="b",pch=19,lwd=2,xaxt="n",col="black",

xlab="month",ylab="value")

axis(1,at=1:12,labels=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))

lines(value[13:24],col="red",type="b",pch=19,lwd=2)

lines(value[25:36],col="orange",type="b",pch=19,lwd=2)

lines(value[37:48],col="purple",type="b",pch=19,lwd=2)

lines(value[49:60],col="blue",type="b",pch=19,lwd=2)

lines(value[61:72],col="green",type="b",pch=19,lwd=2)

legend("bottomright",legend=c("1995","1996","1997","1998","1999","2000"),

lty=1,lwd=2,pch=rep(19,6),col=c("black","red","orange","purple","blue","green"),

ncol=1,bty="n",cex=1.2,

text.col=c("black","red","orange","purple","blue","green"),inset=0.01)

是首先画第一年的，然后逐年添加一条线，然后画图例，算是蛮复杂的。

image003

还有一个简单的方法，就是用ggplot这个包来画。

values <-ozone[1, 1, ]

months=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

years=c("1995","1996","1997","1998","1999","2000")

dat=data.frame(month=rep(months,6),value=values,year=rep(years,each=12))

ggplot(dat,aes(x=month,y=value,group=year,colour=year))+geom_line()

image005

四：测试一下稳健回归模型

稳健回归是加权最小二乘回归，或称文艺最小二乘回归。

MASS 包中的 rlm命令提供了不同形式的稳健回归拟合方式。

回归分析就是用数理统计的方法，研究自然界中变量之间存在的非确定的相互依赖和制约关系，并把这种非确定的相互依赖和制约关系用数学表达式表达出来。其目的在于利用这些数学表达式以及对这些表达式的精度估计，对未知变量作出预测或检验其变化，为决策服务。

介绍几个线性回归（linear regression）中的术语：

残差（Residual）: 基于回归方程的预测值与观测值的差。

离群点（Outlier）: 线性回归（linear regression）中的离群点是指对应残差较大的观测值。也就是说，当某个观测值与基于回归方程的预测值相差较大时，该观测值即可视为离群点。 离群点的出现一般是因为样本自身较为特殊或者数据录入错误导致的，当然也可能是其他问题。

杠杆率（Leverage）: 当某个观测值所对应的预测值为极端值时，该观测值称为高杠杆率点。杠杆率衡量的是独立变量对自身均值的偏异程度。高杠杆率的观测值对于回归方程的参数有重大影响。

影响力点：（Influence）: 若某观测值的剔除与否，对回归方程的系数估计有显著相应，则该观测值是具有影响力的，称为影响力点。影响力是高杠杆率和离群情况引起的。

Cook距离（Cook's distance）: 综合了杠杆率信息和残差信息的统计量。

values <-ozone[1, 1, ]

month.abbr=c("Jan", "Feb", "Mar", "Apr", "May",

"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

month_72 <-factor(rep(month.abbr, length = 72), levels = month.abbr)

deseas1 <-rlm(value ~ month_72 -1)

summary(deseas1)

image007

五．对我们的24*24个地理位置的数据都做稳健回归分析

deseasf<-function(value) rlm(value ~ month_72 -1, maxit= 50)

models <-alply(ozone, 1:2, deseasf) #model是一个list，储存着24*24次稳健回归的结果

failed <-laply(models, function(x) !x$converged)

coefs<-laply(models, coef)#coefs是一个三维数组，记录了所有24×24个位置中每个位置的12个系数

dimnames(coefs)[[3]] <-month.abbr

names(dimnames(coefs))[3] <-"month"

deseas<-laply(models, resid) #deseas是一个三维数组，记录了所有24×24个位置中每个位置的72个残差

dimnames(deseas)[[3]] <-1:72

names(dimnames(deseas))[3] <-"time"

六．对我们的稳健回归系数的三维矩阵进行降维处理，方便画图

通过reshape包可以对三维数组进行降维

coefs_df<-melt(coefs)

head(coefs_df)

lat   long month   value

1 -21.2 -113.8   Jan 264.3964

2 -18.7 -113.8   Jan 261.3284

3 -16.2 -113.8   Jan 260.9643

4 -13.7 -113.8   Jan 258.9999

5 -11.2 -113.8   Jan 255.9999

6 -8.7 -113.8   Jan 254.9999

可以看到第三维的month成功被降维了

还可以通过plyr这个数据工厂包来进行降维

coefs_df<-ddply(coefs_df, .(lat, long), transform, avg= mean(value), std= value / max(value))

>head(coefs_df)   lat   long month   value     avg       std1 -21.2 -113.8   Jan 264.3964 268.6604 0.92770682 -21.2 -113.8   Feb 259.2036 268.6604 0.90948623 -21.2 -113.8   Mar 255.0000 268.6604 0.89473684 -21.2 -113.8   Apr 252.0052 268.6604 0.88422885 -21.2 -113.8   May 258.5089 268.6604 0.90704866 -21.2 -113.8   Jun 265.3387 268.6604 0.9310129

可以看到，不仅成功降维了，还添加了几个属性变量

 

七．最后对降维的coef系数数据画热图

 

coef_limits<-range(coefs_df$value)

coef_mid<-mean(coefs_df$value)

monthsurface<-function(mon)

{

df<-subset(coefs_df, month == mon)

qplot(long, lat, data = df, fill = value, geom="tile") +

scale_fill_gradient(limits = coef_limits,

low = "lightskyblue", high = "yellow")

}

pdf("ozone-animation.pdf", width = 8, height = 8)

l_ply(month.abbr, monthsurface, .print = TRUE)

dev.off()

会在当前R的工作目录下面看到一个pdf文件，里面储存着12个月的在24*24个地理位置的系数热图。

image009