Category Archives: 计算机基础

七 30

用perl把含有简并碱基的引物序列还原成多条序列

Posted on by

这篇博客的程序是错的，请看我最新博客：http://www.bio-info-trainee.com/1528.html

简并碱基对应表格如下；

R:ag
Y:CT
M:AC
K:GT
S:gc
W:AT
H:atc
B:gtc
V:gac
D:GAT
N:ATgc
把这个文本拷贝到txt文件里面保存好，或者直接写入hash里面也行！

[perl]

open FH,"ATCG.txt";

while(<FH>){

chomp;

@F=split/:/;

$hash{$F[0]}=uc $F[1];#右边就是简并表格

}

open FH1,"primer.txt";

while(<FH1>){

next if />/;

chomp;

primer2multiple($_); #对每个含有简并碱基的引物都进行以下函数处理

}

sub primer2multiple{

$primer=$_[0];

$prod=1;

$primer_len=length $primer ;

foreach $i (0..$primer_len-1){

$char=substr($primer,$i,1);

if ($char !~/[ATCG]/){$prod*=length $hash{$char}}

}

$new="";

foreach $i (0..$primer_len-1){

$char=substr($primer,$i,1);

if ($char =~/[ATCG]/){$new.=$char x $prod}

else {$tmp=length $hash{$char};$new.=$hash{$char} x ($prod/$tmp)}

}

die "error!" if   $primer_len*$prod != length $new ;

foreach $i (0..$prod-1){

$tmp="";

for(my $j=$i;$j<(length($new));$j+=$prod){$tmp.=substr($new,$j,1)}

print "$tmp\n";

}

}

[/perl]

可以直接调用这个函数即可primer2multiple('ATGCVHGT');

就可以看到简并碱基被替换啦，就是那个V和H

 

ATGCGAGT
ATGCATGT
ATGCCCGT
ATGCGAGT
ATGCATGT
ATGCCCGT
ATGCGAGT
ATGCATGT
ATGCCCGT

七 30

perl用递归法得到允许多个错配的模式匹配

Posted on by

如果我要匹配一个字符串$a="ATTCCGGGAT";那么直接在shell里面grep它即可，写脚本也行$seq =~ /$a/;，但是如果我想查找这个字符串的模糊匹配，允许一个错配的情况，那么就非常多了！这时候简单的匹配已经不能达到目的，但是我们仍然可以用perl强大的正则匹配功能达到目的。

比如，它的匹配模式应该:

/.TTCCGGGAT/

/A.TCCGGGAT/

/AT.CCGGGAT/等等，可以把这些模式都综合起来就是下面这个

$b=(?^:(?:A(?:T(?:T(?:C(?:C(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT)|.GGGAT)|.CGGGAT)|.CCGGGAT)|.TCCGGGAT)|.TTCCGGGAT))

所以我们就应该通过程序来生成这个字符串，然后用

$seq =~ /$b/;来替代$seq =~ /$a/;

而允许两个错配的格式就更复杂了：

(?^:(?:A(?:T(?:T(?:C(?:C(?:G(?:G(?:G..|.(?:A.|.T))|.(?:G(?:A.|.T)|.AT))|.(?:G(?:G(?:A.|.T)|.AT)|.GAT))|.(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT))|.(?:C(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT)|.GGGAT))|.(?:C(?:C(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT)|.GGGAT)|.CGGGAT))|.(?:T(?:C(?:C(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT)|.GGGAT)|.CGGGAT)|.CCGGGAT))|.(?:T(?:T(?:C(?:C(?:G(?:G(?:G(?:A.|.T)|.AT)|.GAT)|.GGAT)|.GGGAT)|.CGGGAT)|.CCGGGAT)|.TCCGGGAT)))

[perl]

$a="ATTCCGGGAT";
$one_match=fuzzy_pattern($a,1);
print "$one_match\n";

sub fuzzy_pattern {
my ($original_pattern, $mismatches_allowed) = @_;
$mismatches_allowed >= 0
or die "Number of mismatches must be greater than or equal to zero\n";
my $new_pattern = make_approximate($original_pattern, $mismatches_allowed);
return qr/$new_pattern/;
}
sub make_approximate {
my ($pattern, $mismatches_allowed) = @_;
if ($mismatches_allowed == 0) { return $pattern }
elsif (length($pattern) <= $mismatches_allowed)
{ $pattern =~ tr/ACTG/./; return $pattern }
else {
my ($first, $rest) = $pattern =~ /^(.)(.*)/;
my $after_match = make_approximate($rest, $mismatches_allowed);
if ($first =~ /[ACGT]/) {
my $after_miss = make_approximate($rest, $mismatches_allowed-1);
return "(?:$first$after_match|.$after_miss)";
}
else { return "$first$after_match" }
}
}

[/perl]

只需要控制$one_match=fuzzy_pattern($a,1);里面的参数即可控制自己想要的匹配情况。

然后把生成的匹配模式用了进行序列匹配$seq =~ /$one_mismatch/;

这个程序的重点就是解析需要生成的匹配字符串规则，然后用递归来生成这个匹配字符串。

这种匹配，在引物搜索特别有用。

七 24

perl多重循环生成全排列精简

Posted on by

我想输出ATCG四个字符，组成一个12个字符长度字符串的全排列。共4^12=16777216种排列，

AAAAAAAAAAAA

AAAAAAAAAAAT

AAAAAAAAAAAC

AAAAAAAAAAAG

AAAAAAAAAATA

`````````````````````````

按照正常的想法是通过多重循环来生成全排列，所以我写下了下面这样的代码，但是有个问题，它不支持多扩展性，如果100个全排列，那么得写一百次循环嵌套。

4

所以我求助了perl-china群里的高手，其中一个给了一个结合二进制的巧妙解决方案。

perl -le '{$a{"00"}="A";$a{"01"}="T";$a{"10"}="C";$a{"11"}="G";for(0..100){print join"",map{$a{$_}} unpack("a2"x12,sprintf("%024b\n",$_))}}'

AAAAAAAAAAAA

AAAAAAAAAAAT

AAAAAAAAAAAC

AAAAAAAAAAAG

AAAAAAAAAATA

AAAAAAAAAATT

AAAAAAAAAATC

AAAAAAAAAATG

AAAAAAAAAACA

AAAAAAAAAACT

AAAAAAAAAACC

AAAAAAAAAACG

AAAAAAAAAAGA

AAAAAAAAAAGT

AAAAAAAAAAGC

AAAAAAAAAAGG

AAAAAAAAATAA

AAAAAAAAATAT

这里我简单解释一下这个代码

perl -le '{$a=sprintf("%024b\n",1);print $a}'

000000000000000000000001

这个sprintf是根据024b的格式把我们的十进制数字转为二级制，但是补全为24位。

perl -le '{@a=unpack("a2"x12,sprintf("%024b\n",100));print join":", @a}'

00:00:00:00:00:00:00:00:01:10:01:00

unpack这个函数很简单，就是把二进制的数值拆分成字符串数组，两个数字组成一个字符串。

最后通过map把$a{"00"}="A";$a{"01"}="T";$a{"10"}="C";$a{"11"}="G";对应成hash值并且输出即可！

然后这个大神又提成了一个递归的办法来解决

[perl]

#!/usr/bin/perl

my $arr = [];

$arr->[$_] = 0 for (0..11);
my @b = ("A","C","G","T");

while(1){
print join "", map {$b[$_]} @$arr;
print "\n";
exit unless &count($arr,11);
}

sub count {
my $arr = shift;
my $x = shift;

return 0 if $x < 0;

if ($arr->[$x] < 3){
$arr->[$x]++;
return 1;
}

else{
$arr->[$x] = 0;
return &count($arr,$x - 1);
}
}

[/perl]

还有另外一个大神也提成了一个递归的算法解决，算法之道还是蛮有趣的

2

最后还补充一个也是计算机技巧的解决方案，也是群里的朋友提出来的。

[perl]

#!/usr/bin/perl -w

my @A = qw(A T G C);
my $N = 12;

foreach my $n (0..4**$N){
foreach my $index (0..$N){
print $A[$n%4];
$n = $n >>2;
}
print "\n";
}

[/perl]

这个是位运算，如果C语言基础学的好的同学很快就能看懂的。

 

七 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

perl操作pdf文档

Posted on by

大家看看就好，这个模块写的不怎么样，而且有高手已经写了一个pdftoolkit就是完全用这个模块实现了大部分pdf文档的操作

PDF::API2模块使用笔记

一：简单使用方法

use PDF::API2;

# Create a blank PDF file $pdf = PDF::API2->new();
# Open an existing PDF file $pdf = PDF::API2->open('some.pdf');
# Add a blank page $page = $pdf->page();
# Retrieve an existing page $page = $pdf->openpage($page_number);
# Set the page size $page->mediabox('Letter');
# Add a built-in font to the PDF $font = $pdf->corefont('Helvetica-Bold');
# Add an external TTF font to the PDF $font = $pdf->ttfont('/path/to/font.ttf');
 

# Add some text to the page

 $text = $page->text();

$text->font($font, 20);

$text->translate(200, 700);

$text->text('Hello World!');
# Save the PDF $pdf->saveas('/path/to/new.pdf');

 

实例：

use PDF::API2;

$pdf=PDF::API2->new;

$pdf->mediabox('A4');

$ft=$pdf->cjkfont('Song');

$page = $pdf->page;

$gfx=$page->gfx;

$gfx->textlabel(50,750,$ft,20,"\x{Cool44}\x{4EA7}"); # 资产二字

$pdf->saveas('Song_Test.pdf');

 

二：主要对象及方法

1、pdf对象可以创造，可以打开，可以保存，可以更新，还有一堆参数可以设置

$pdf->preferences(%options)还可以设置一些浏览参数，不过本来pdf阅读器可以设置，没必要在这里花时间。

这个可以当做是个人创建pdf的保密信息，也许有一点用吧。

还可以可以设置页脚$pdf->pageLabel($index, $options

2、Page对象，可以新建，可以打开，可以保存（需要指定保存的位置）

$page = $pdf->page()

$page = $pdf->page($page_number)

$page = $pdf->openpage($page_number);

还可以更新旧的pdf，这样可以循环获取pdf页面不停的累积到一个新的pdf

$page = $pdf->import_page($source_pdf, $source_page_number, $target_page_number)

$pdf = PDF::API2->new();

$old = PDF::API2->open('our/old.pdf');   # Add page 2 from the old PDF as page 1 of the new PDF

$page = $pdf->import_page($old, 2);

$pdf->saveas('our/new.pdf');If $source_page_number is 0 or -1, it will return the last page in the document.

$count = $pdf->pages()Returns the number of pages in the document.

这样就可以写一个简单程序把我们的pdf文件合并

use PDF::API2;

my $new = PDF::API2->new;

foreach my $filename (@ARGV) {   my $pdf = PDF::API2->open($filename);   $new->importpage($pdf, $_) foreach 1 .. $pdf->pages;}$new->saveas('new.pdf'); $pdf->mediabox($name)

可以指定A4，A3，A5等等$pdf->mediabox($w, $h)可以指定宽度和高度$pdf->mediabox($llx, $lly, $urx, $ury)

3，还可以随意画点线面及表格，太复杂了就不看了

 

七 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

七 15

perl操作excel表格

Posted on by

perl这个语言已经过时很久了，所以它的模块支持性能不是很好，暂时我只看到了对excel2003格式的表格的读取写入操作。

以下是我参考Spreadsheet::ParseExcel这个模块写的一个把excel表格转为csv的小程序，大家也可以自己搜索该模块的说明文档，这样学的更快一点！

[perl]

#!/usr/bin/perl -w
# For each tab (worksheet) in a file (workbook),
# spit out columns separated by ",",
# and rows separated by c/r.

use Spreadsheet::ParseExcel;
use strict;
use utf8;
use Encode::Locale qw($ENCODING_LOCALE_FS);
use Encode;

my $filename ="test.xls";#输入需要解析的excel文件名，必须是03版本的
my $e = new Spreadsheet::ParseExcel; #新建一个excel表格操作器
my $eBook = $e->Parse($filename);    #用表格操作器来解析我们的文件
my $sheets = $eBook->{SheetCount};   #得到该文件中sheet总数
my ($eSheet, $sheetName);

foreach my $sheet (0 .. $sheets - 1) {
$eSheet = $eBook->{Worksheet}[$sheet];
$sheetName = $eSheet->{Name};
my $f1 = encode(locale_fs => $sheetName); #每次操作中文我都很纠结，还得各种转码
open FH_out ,">$f1.csv" or die "error open ";
next unless (exists ($eSheet->{MaxRow}) and (exists ($eSheet->{MaxCol})));
foreach my $row ($eSheet->{MinRow} .. $eSheet->{MaxRow}) {
foreach my $column ($eSheet->{MinCol} .. $eSheet->{MaxCol}) {
if (defined $eSheet->{Cells}[$row][$column])
{
print FH_out $eSheet->{Cells}[$row][$column]->Value . ",";
} else {
print FH_out ",";
}
}
print FH_out "\n";
}
close FH_out;

}
exit;

[/perl]

 

五 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