在一些天文期刊上经常看到赤道坐标星图，不但给出了某一个时刻太阳、月亮的位置，还给出了大行星在一年中相对于星空背景的运行轨迹。今天我们就用R软件绘制这样的星图。

当然，为了绘制星图，我们需要知道恒星的位置和亮度。这就需要引入天文坐标系。天文学上天体的位置有多种坐标系统，包括赤道坐标系，黄道坐标系和地平坐标系等。

首先介绍一下今天用到的赤道坐标系。

赤道坐标系是以春分点作为起点，春分点和秋分点是赤道和黄道两个大圈相交的点。由于地球的自转，我们看到太阳相对于恒星背景从西向东运行，太阳所经过的路线就称为黄道。不难想象，春分点过后，太阳相对于星空背景一天天逐渐向北移动（当然也同时向东移动，差不多每天移动1度多一点儿）。我们就说春分点为升交点。太阳到达夏至点的时候，是最北的，之后逐渐向南，到秋分点又与赤道相交。因此秋分点又称为降交点。太阳继续向南，此时太阳位于天赤道以南。到达冬至点之后又逐渐向北。如此周而复始。

我们定义赤道坐标系的原点就是升交点，向左度量，一圈用24h表示。其经度就称为“赤经”，其纬度就称为“赤纬”，用度表示， 北半球为正，南半球为负。这就是我们将要看到的图形中的星图。天文学上的星图左东右西，与地图正好相反。

其次，我们需要了解一下恒星的位置和亮度。
这里我们选取了耶鲁亮星星表。这个星表包含了肉眼可见的9110颗较亮恒星的位置亮度等信息。按照肉眼可见的亮度分成不同的星等，星的亮度是用照相测量法测量的。

第三，我们需要知道行星在任何一个时刻的位置。即，我们输入时刻，就有程序可以计算出该天体的位置等信息。这是球面天文学和天体力学的内容了。我们这里使用R的moonsun程序包计算某一时刻天体的位置和亮度。 我们需要输入的时间为儒略日（Julian day）。儒略日是什么？儒略日是指从公元前4713年1月1日中午12时开始，开始距离现在的天数。我们输入的日期都必须先转换成儒略日，再用儒略日作为时间点计算。儒略日计算两个时间点之间相隔的天数非常方便，因此在天文学上得到了广泛的应用。

第四，我们需要了解一下行星运行的基本知识。
水星，金星位于地球内侧，更靠近太阳，称为内行星。火星，木星，土星，天王星，海王星位于地球外侧，称为外行星。
八大行星的轨道与黄道相差不大，而且均自西向东运行。靠近太阳最近的水星运行速度最快。从地球上看，行星相对于星空背景有时候向东，有时候向西。向东称为顺行，各行星大部分时间在顺行中，向西称为逆行，在顺行和逆行的转换过程中，有较短的时间，行星相对星空背景处于相对静止的状态，称为留。

这里给出了skyatlasR程序能够绘制行星周年视运动的轨迹，当然，由于没有考虑到星体之间的摄动等更精确的天体力学因素，行星的位置计算不太高，精度在几角分之内，但是仍然能够满足大部分天文科普的需求。更精确的历表就得参考NASA,也就是美国国家航空航天局喷气推进实验室的HORIZON历书系统了。用户可以在线查询一些信息。

下面简要介绍一下 skyatlas 程序的用法

使用方法

当然，首先要下载我们编好的程序了。 下载 skyatlas.rar

1 下载R软件 http://ftp.ctex.org/mirrors/CRAN/bin/windows/base/R-2.11.1-win32.exe 下载
2 将 “plot.r”, “resstar.csv”, “moonsun_0.1.1.zip” 拷贝到 C:下,”moonsun_0.1.1.zip” 为R程序包，安装时不需要解压缩。
3 安装moonsun.zip程序包，双击R的图标，路径为 开始>所有程序> R, 2.11.1, 然后再菜单中选择，菜单“程序包”，“从本地zip安装程序包”，选择 moonsun_0.1.1.zip 即可。
4 从R的菜单 “文件” “source R code”. 打开 “C:/plot.r”,或者输入命令source(“C:/plot.r”)即可。

运行实例
R图形的组件都是一步一步添加上去的，故绘制一个图形步骤较多，这里举出了绘制水星1980年和2011年的例子，可供有兴趣的读者参考。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

######################################
######################################

### 例一 绘制1980年水星的周年视运动  
##  该结果与《中国大百科全书 天文学卷》参见 485页 进行比较，水星的运行曲线结果是一致的。
x11(14, 5)
##绘制底图
background(type = "n", xlim =c(0, 24), ylim = c(-40, 40), xlab = "赤经", ylab = "赤纬")

##添加黄道
add.ecliptic()

##添加4等以上的恒星
add.stars(resstar, mag = 4, pch = 19, time = 1)

## 添加水星的运行轨迹
tracks(mercury(jd(year = 1980,month = 1,day = 1,length = 365)))
## 添加亮星的名称
add.bright.lab ()

###生成每隔30天的点数据，
position <- mercury(jd(year = 1980,month = 1, day = 1,length = 12, by = 30))
addlab(position)
title("1980年水星的周年视运动")

结果如下

例二

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

### 2011水星的周年视运动

### 第一步 打开一个绘图窗口
x11(14, 5)
##  第二步 绘制底图
background(type = "n", xlim =c(0, 24), ylim = c(-40, 40), xlab = "赤经", ylab = "赤纬")

##添加黄道
add.ecliptic()

##第四步 添加4等以上的恒星
add.stars(resstar, mag = 4, pch = 19, time = 1)

## 第五步 添加水星的运行轨迹
tracks(mercury(jd(year = 2010,month = 1,day = 1,length = 365)))

## 第六步 添加亮星的名称
add.bright.lab()

###第七步 添加生成每隔30天的点数据，并标在图上
position <- mercury(jd(year = 2010,month = 1, day = 1,length = 12, by = 30))
addlab(position)

### 第八步 添加标题
title("2011年水星的周年视运动")

例三 当前时刻日月及大行星的位置图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

### 打开一个绘图窗口
x11(14, 5)

## 绘制赤道背景坐标
background(type = "n", xlim =c(0, 24), ylim = c(-40, 40), xlab = "赤经", ylab = "赤纬", bg = "darkblue")

##添加黄道
add.ecliptic()

## 添加4等以上的恒星
add.stars(resstar, mag = 4, pch = 19, time = 1)

## 添加亮星的名称
add.bright.lab()

#### 添加大行星
add.planet(col = "black")

## 添加标题
title(paste("", "日月及大行星的位置(",as.character(Sys.Date()),")" ))

生态位理论是解释生物多样性格局和物种共存的重要理论之一。每个物种对生境的偏好，就确定了其生态位。从1910年至1980年以前，国外生态学家对物种生态位的计算等生态位的宽度和计算进行了大量研究，主要的计算方法都是那时候提出的。从上世纪80年代到现在，我国的生态学界仍然能见到一些计算物种生态位宽度和生态位重叠的论文。

张金屯老师的《数量生态学》更是对生态位和生态位重叠的计量方法进行了详细介绍。

为了方便的完成这些指数的计算，本人粗略编写了一些R程序，完成相应指数的计算.还增加了Bootstrap检验，从而获得生态位重叠的置信区间估计。

输入数据的格式是标准的群落物种矩阵格式。即每一行代表每一个样方，每一列代表每一个物种。

下面就是相应的R程序，可供下载。

由于当前代买还存在很多冗余，欢迎各位同行提出宝贵意见。以便修订之后整合到进行种间关联度分析的spaa 程序包中。

有任何意见或建议，均可发邮件至 jinlongzhang01@gmail.com 来信必复。

种间关联度及生态位重叠的R程序包

注意！ 运行Lagrange软件的指南随后贴出，现在的版本仅供参考。

“Lagrange is a Python package implementing likelihood models for geographic range evolution on phylogenetic trees, with methods for inferring rates of dispersal and local extinction and ancestral ranges. “ http://code.google.com/p/lagrange/

1 Download and install Python NumPy and SciPy

Lagrange requires Python, NumPy and SciPy, the first step is to download and install these software.

Download and install Python2.5, because SciPy requires this version, it can be downloaded fromhttp://www.python.org/download/

Download and install NumPy and SciPy ,you can down load the exe files form

http://sourceforge.net/project/showfiles.php?group_id=27747 for SciPy

http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=175103 for NumPi

It is strongly recommended that you download the numpy-1.2.1-win32-superpack-python2.5.exe, since the imcompatible problems may exist if you installed a incorrect version.

You can visit scipy’s website for more detail.

http://www.scipy.org/

2 Copy the Lagrange modules to working directory

Assume you have had all the software that Lagrange required installed on your PC. You may seed a directory under C:Python25

now unzip your lagrange2[1].0.1.tar.gz

copy all the file in the directory which just unzipped to C:python25

3 Data preparation

The python scripts for exercise can be downloaded at

http://bodegaphylo.wikispot.org/Biogeography_%28Donoghue_%26_Smith%29

copy the test python script named psychotria.py and psychotria.matrix to python directory, C:\python25 .

4 Run python scripts in cmd

The python scripts can be implement in commander just which just like the DOS , to run the script, begin here:

Begin>run> type cmd

a commander window will appeared

just type cd C:python25 to change your working directory to test

and type psychotria.py>results.txt “>results.txt” allowed you to save the results of the scripts in “results.txt”

When all are got ready , the python script will be run in your PC

1. seqRFLP简介

RFLP和TRFLP等技术在生物多样性研究中十分广泛。随着测序技术的日趋普及，在进行RFLP或TRFLP之前，用户很可能已经获得了DNA序列。为了对RFLP和TRFLP的结果进行更好的分析，特别是选择合适的内切酶，对模糊的电泳条带进行确认，基于RFLP/TRFLP及序列分别对物种进行鉴定等显得日趋重要，但是目前还缺乏这样一种已知序列就能获得RFLP/TRFLP的分析平台。为此我们编写了seqRFLP软件。

seqRFLP是用来对DNA序列进行模拟酶切、模拟PCR引物筛选的软件包，目前版本是1.0.0。用开源的R语言 (http://www.r-project.org/ )写成，可以在Linux, Windows， MacOS等多种计算机平台上运行。seqRFLP已经被RCRAN接收，用户可以在R的任何一个镜像下载。本说明假设用户使用Windows，其他平台用户的使用方法与之类似，均需先安装R软件。

2. seqRFLP的下载和安装

首先下载和安装R软件。建议选择较新的版本，如Windows: R-2.11.1，其下载地址之一为:

http://ftp.ctex.org/mirrors/CRAN/bin/windows/base/R-2.11.1-win32.exe

下载到本地后，双击R-2.11.1-win32.exe开始安装，各选项默认即可。安装完成后双击

开始>所有程序>R>R 2.11.1

输入命令:

1

install.packages("seqRFLP")

R将提示选择CRAN,选择距离较近的CRAN镜像，将软件包自动下载和安装好。

导入seqRFLP程序包：

1

library(seqRFLP)

3.查询seqRFLP的帮助

seqRFLP 的每个函数都有详细的帮助文件，要查询帮助可输入:

?seqRFLP 查看详尽的帮助手册:用户也可以在CRAN上下载pdf版本的帮助手册。

seqRFLP符合GPL-2协议，用户可以免费使用并拷贝该软件，也可以更改其源代码。

运行seqRFLP实例:

1

example(seqRFLP)

致 谢

感谢马克平研究员、梁宇博士的指导，感谢博士研究生孙秀峰及龙恩熙对软件的测试，感谢软件包的编写有过帮助的研究组各同学和老师。

参考文献

Saiki RK, Scharf S, Faloona F, Mullis KB, Erlich HA, Arnheim N (1985). Enzymatic amplification of beta-globin genomic sequences and restriction site analysis for diagnosis of sickle cell anemia. Science 230 (4723):1350-4

Roberts, R.J., Vincze, T., Posfai, J., Macelis, D. (2010) REBASE–a database for DNA restriction and modification: enzymes, genes and genomes. Nucl. Acids Res. 38: D234-D236. http://rebase.neb.com

Phylomatic (http://phylodiversity.net/phylomatic/) 和PhyloCom (http://phylodiversity.net/phylocom/) 是Cam Webb、Steve Kembel和David Ackerly编写的软件，用于群落系统发育分析。

Phylomatic可基于植物名录以及相应进化树骨架生成进化树。不过，通过Phylomatic建立的进化树，一般只能准确反映科或者属以上分类等级的进化关系（取决于该属是否在骨架进化树上是否有信息）。由于不能准确反映科或属以下的系统发育关系，Phylomatic生成的进化树，一般只用于群落系统发育，或者空间尺度较大，类群较多，且系统发育关系较远的系统发育多样性分析中。如果分类单元数比较少，在Genbank中又有足够的序列， 则应该优先考虑自行建立分子系统树，再进行相应分析，具体请参考 （Roquet et al. 2013) 。

通过植物名录建立进化树， 目前也有两个程序包可以选择， 一个是浙江大学金毅老师编写的 S.PhyloMaker R 脚本 (https://github.com/jinyizju/S.PhyloMaker), 另一个是 Scott Chamberlain 编写的 brranching R 程序包 (https://github.com/ropensci/brranching)。 S.PhyloMaker 的数据是基于Zanne 2014进化树生成的，可建立包含几万个分类单元的进化树，但速度较慢。 brranching 主要是基于Phylomatic的，可以直接输入植物名录获得进化树，速度较快，但是在分类单元数量上有一定限制。以上两个程序包的用法，用户可自行参考指南，本文只介绍通过Phylomatic建立进化树。

Phylocom是用来进行群落系统发育分析的软件。该软件可以为Phylomatic软件得到的无枝长进化树，用BLADJ方法拟合出枝长（这是因为Phylomatic软件建立的进化树， 有科属种的标签，而另外一个提供的ages文件， 有科属种的分化时间），不过BLADJ算法在2014年以后已经很少使用。此外，Phylocom还可以：计算系统发育多样性（PD），系统发育结构(community structure)，群落的系统发育距离(community phylogenetic distance)，分析群落的性状进化(AOT)等。下面就介绍如何使用Phylomatic和Phylocom进行相应的分析。

假设有如下植物名录要建立进化树：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

Castanopsis eyrei
Castanopsis carlesii
Castanopsis fargesii
Castanopsis tibetana
Cyclobalanopsis gracilis
Cyclobalanopsis glauca
Lithocarpus glaber
Myrica rubra
Elaeocarpus decipiens
Syzygium buxifolium
Daphniphyllum oldhamii
Loropetalum chinense
Rhododendron latoucheae
Rhododendron ovatum
Vaccinium carlesii
Schima superba
Adinandra millettii
Eurya muricata
Camellia fraterna
Ternstroemia gymnanthera
Machilus thunbergii
Neolitsea aurata var. chekiangensis
Chimonanthus salicifolius

一般流程如下

1. 准备 科/属/种名录

请注意: 《中国植物志》 以及Flora of China等植物志, 由于成书较早，所以采用的分类系统都不是APG。而APG分类系统是目前为止能够较为客观体现科水平系统发育关系的分类系统（最新版本APG IV），因而是推荐采用的。不过目前能获得的准确科属列表为APGIII的版本，与APG IV差别不大。为了批量查询植物所在APG的科名，可使用本人编写的R程序包 plantlist (https://github.com/helixcn/plantlist)。该程序包的简要介绍参见 http://blog.sciencenet.cn/blog-255662-846673.html 。

注意： 由于plantlist给出的是APGIII的科属列表，如果用户要使用APG IV，请自行修订。

1.1 plantlist程序包的安装

1
2

library(devtools)
install_github("helixcn/plantlist")

1.2 生成科属种列表

R 命令为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

sp <- c( "Castanopsis eyrei","Castanopsis carlesii",
"Castanopsis fargesii","Castanopsis tibetana",
"Cyclobalanopsis gracilis","Cyclobalanopsis glauca",
"Lithocarpus glaber","Myrica rubra",
"Elaeocarpus decipiens","Syzygium buxifolium",
"Daphniphyllum oldhamii","Loropetalum chinense",
"Rhododendron latoucheae","Rhododendron ovatum",
"Vaccinium carlesii","Schima superba",
"Adinandra millettii","Eurya muricata",
"Camellia fraterna","Ternstroemia gymnanthera",
"Machilus thunbergii","Neolitsea aurata var. chekiangensis",
"Chimonanthus salicifolius")

library(plantlist)
sp2 <- TPL(sp)
res <- taxa.table(sp2)
writeLines(res, "taxa.table.txt")
getwd()

运行以上R代码，科属种列表就保存在taxa.table.txt 文件中，保存的路径在R中有所显示。

所生成文件taxa.table.txt用Notepad++ (https://notepad-plus-plus.org/download/ ) 记事本打开， 内容如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

Fagaceae/Castanopsis/Castanopsis_eyrei
Fagaceae/Castanopsis/Castanopsis_carlesii
Fagaceae/Castanopsis/Castanopsis_fargesii
Fagaceae/Castanopsis/Castanopsis_tibetana
Fagaceae/Cyclobalanopsis/Cyclobalanopsis_gracilis
Fagaceae/Cyclobalanopsis/Cyclobalanopsis_glauca
Fagaceae/Lithocarpus/Lithocarpus_glaber
Myricaceae/Myrica/Myrica_rubra
Elaeocarpaceae/Elaeocarpus/Elaeocarpus_decipiens
Myrtaceae/Syzygium/Syzygium_buxifolium
Daphniphyllaceae/Daphniphyllum/Daphniphyllum_oldhamii
Hamamelidaceae/Loropetalum/Loropetalum_chinense
Ericaceae/Rhododendron/Rhododendron_latoucheae
Ericaceae/Rhododendron/Rhododendron_ovatum
Ericaceae/Vaccinium/Vaccinium_carlesii
Theaceae/Schima/Schima_superba
Theaceae/Camellia/Camellia_fraterna
Pentaphylacaceae/Adinandra/Adinandra_millettii
Pentaphylacaceae/Eurya/Eurya_muricata
Pentaphylacaceae/Ternstroemia/Ternstroemia_gymnanthera
Lauraceae/Machilus/Machilus_thunbergii
Lauraceae/Neolitsea/Neolitsea_aurata_var._chekiangensis
Calycanthaceae/Chimonanthus/Chimonanthus_salicifolius

注： 查询植物科名，也可以用 taxize程序包（https://ropensci.org/tutorials/taxize_tutorial/）

2 Phylomatic建树

在Zanne 2014进化树发表之前，用Phylomatic建立进化树一般只是按照APGIII的骨架，生成无枝长的newick进化树。然后再在Phylocom软件中，用BLADJ模块，将Wikstroem等人推断的各类群分化时间拟合到无枝长的进化树上。Wikstroem的分化时间由于发表年代久远， 被子植物各类群的分化时间有一些变动，通过Phylocom BLADJ的方法校正分化时间已经显得过时。从2015年开始, Phylomatic网站整合了 Zanne 2014的进化树骨架, 可以直接基于该进化树直接获得有枝长的进化树。

方法如下:

1. 打开 http://phylodiversity.net/phylomatic/
2. 将 科/属/种 列表 拷贝到 taxa = 框中
3. storedtree = 选择 Zanne 2014
4. 勾选 clean = TRUE， 其余选择默认.
5. 点击Send按钮获得结果

将网页中的输出结果拷贝到文本文件中， 并另存为 phylo文件， 请注意，如果后续用Phylocom进行分析，则该phylo文件不能有任何扩展名。该进化树为Newick格式， 不含Singletons，因此可以用R的ape读取并进行数据处理，或者Figtree软件绘图。

获得的进化树如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

((Chimonanthus_salicifolius:108.8,(Neolitsea_aurata_var._chekiangensis:88.1284,
Machilus_thunbergii:88.1284)Lauraceae:20.672)Laurales:71.8279,
(((((((Cyclobalanopsis_gracilis:20.205,Cyclobalanopsis_glauca:20.205)
Cyclobalanopsis:20.205,(Lithocarpus_glaber:3.30372,
(Castanopsis_eyrei:0.94544,(Castanopsis_carlesii:0.320126,
(Castanopsis_tibetana:0.28527,Castanopsis_fargesii:0.28527):0.034856):
0.625313):2.35828):37.1063)Fagaceae:21.8954,Myrica_rubra:62.3054):
50.9532,Elaeocarpus_decipiens:113.259)Fabidae:4.27613,Syzygium_buxifolium:
117.535)Rosidae:0.138911,(Daphniphyllum_oldhamii:107.877,Loropetalum_chinense:
107.877):9.79685)Superrosidae:1.49186,((Rhododendron_latoucheae:62.0517,
Vaccinium_carlesii:62.0517,Rhododendron_ovatum:62.0517)Ericaceae:20.7165,
((Ternstroemia_gymnanthera:43.8733,(Adinandra_millettii:33.7061,
Eurya_muricata:33.7061):10.1672):27.9519,(Schima_superba:57.076,
Camellia_fraterna:57.076):14.7492):10.943):36.3974)
Gunneridae:61.4628):259.665;

注：用brranching R程序包可以直接输入物种名，获得进化树。

3. 用Phylocom计算PD, NRI 和NTI （以Example文件为例）

注意：这部分计算使用的是Example文件，与上文的植物名录已经没有关系。

3.1 下载phylocom软件

打开 https://github.com/phylocom/phylocom/releases 下载zip文件，解压缩, 可见以下几个文件夹：

1. Example_data 中是练习数据
2. Mac 中是苹果机的运行软件
3. R，是驱动Phylocom的R脚本
4. Src 是Phylocom的源程序，phylocom是C语言写的
5. W32 是Windows平台下可以运行的exe程序
6. Phylocom_manual.pdf 是软件说明书
7. README 注意事项

3.2 创建工作路径

Phylocom必须在纯ASCII码的路径下运行，运行路径不能有中文。在C盘根目录下，创建一个名为phylocom的文件夹。将w32中的phylocom.exe文件，拷贝到文件夹 C:/phylocom下。 将example_data文件夹中的phylo和sample文件，拷贝到该文件夹C:/phylocom中。

3.3 利用phylocom软件计算PD、NRI、NTI等指数

sample文件第一列是样方的名称，第二列是个体数，第三列是物种名。物种名一定要与phylo文件中的物种名完全对应（包括大小写）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

clump1    1    sp1
clump1    1    sp2
clump1    1    sp3
clump1    1    sp4
clump1    1    sp5
clump1    1    sp6
clump1    1    sp7
clump1    1    sp8
clump2a    1    sp1
clump2a    1    sp2
.... (中间省略)
random    2    sp15
random    3    sp17
random    1    sp22
random    2    sp24

将phylo和sample文件夹拷贝到 C:\phylocom 之后，

点击 开始> cmd > 输入 cd C:\phylocom>

输入 phylocom comstruct > comstruct.output 计算群落的系统发育结构的NRI和NTI指数

输入 phylocom pd > pd.output.txt 计算系统发育多样性 Faith’s PD

输入 phylocom comdist > comdist.output.txt 计算群落两两之间系统发育距离 MPD的MNDT等

注意： 进行系统发育多样性相关的计算， 可以参考picante 程序包 (http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=1116212)

参考文献

1. Chamberlain, S. (2018). brranching: Fetch ‘Phylogenies’ from Many Sources. R package version 0.3.0. https://CRAN.R-project.org/package=brranching
2. Qian, H. and Y. Jin. (2016) An updated megaphylogeny of plants, a tool for generating plant phylogenies and an analysis of phylogenetic community structure. Journal of Plant Ecology 9(2): 233–239.
3. Roquet, C., Thuiller, W., & Lavergne, S. (2013). Building megaphylogenies for macroecology: taking up the challenge. Ecography, 36(1), 13-26.
4. Webb, C. O., & Donoghue, M. J. (2005). Phylomatic: tree assembly for applied phylogenetics. Molecular Ecology Notes, 5(1), 181-183.
5. Zhang, J.L. (2018). plantlist: Looking Up the Status of Plant Scientific Names based on The Plant List Database. R package version 0.3.7/r31. https://R-Forge.R-project.org/projects/plantlist/

编写一个计算种间关联的小程序是出于以下原因：

查阅了很多进行种间关联分析的论文，一直也没有提及用什么软件。 为此我觉得有必要编写一个程序包，用来进行群落数据的处理及种间关联度计算和绘图。

R当然成为首选，因为R是开源的，所有代码可见；同时又是免费的，不用注册，就可以使用R。

花费几天的时间，终于完成了。定名为spaa, 意为 species association analysis。可以在R的CRAN上下载

http://cran.r-project.org/web/packages/spaa/index.html

spaa中有如下几个函数：

• data2mat 将物种列表转换为矩阵
• freq.calc 计算物种频度
• plot.lowertri 绘制半矩阵图
• plot.network 绘制种间关联图 （有待进一步完善）
• sp.assoc 计算种间总体关联性
• sp.pair 计算物种两两之间的关联性
• sub.sp.matrix 在总矩阵中，按照物种出现的频度筛选物种。

群落按照物种相似形组成进行聚类分析，可以用树状图较好的表现物种的组成关系。受到很多植被学家的重视。这里以R软件实现聚类分析为例。

如果按照物种组成的相似性做聚类分析，那么可以用Jaccard指数（经过转换的）。Jaccard指数只考虑物种在两个样方间是否重复出现，盖度在分析的过程中并不起什么作用。但是如果对乔木和灌木进行分析，就可以考虑个体的数量，计算样方物种组成的相似性的时候用Bray-Curtis指数。Jaccard指数和Bray-Curtis指数在众多生态学相关的程序包中都是可以计算的。下面说一下在R软件中，结合vegan程序包，对草本样方的物种组成进行聚类分析。

下面是在R中的具体操作过程：

第一步 矩阵的整理，建议先整理一下各样地的名录，成如下格式，再用R整理成物种矩阵。

1
2
3
4
5
6
7
8
9
10
11
12
13

plotname species
plot1 sp1
plot1 sp2
plot1 sp3
plot1 sp4
plot1 sp5
plot2 sp1
plot2 sp3
plot3 sp4
plot3 sp2
plot3 sp6
plot3 sp7
.....

在Excel中，另存为csv格式，如存名称为 herbplots.csv。

第二步 读取文件

1

herb.data<- read.csv("D:/herb/herbplots.csv", header=T)

第三步 转换为 矩阵

1

pre.matrix <- table(herb.data)

此时生成的矩阵，形式如下：

1
2
3
4

plotname sp1 sp2 sp3 sp4 sp5 sp6 sp7
plot1 1 1 1 1 1 0 0
plot2 1 0 1 0 0 0 0
plot3 0 1 0 1 0 1 1

第四步 计算各样方的Jaccard距离

1

plot.dist <- vegdist(pre.matrix, method="jaccard")

距离矩阵为一个半矩阵，形式如下

1
2
3

plot1 plot2
plot2 0.6000000
plot3 0.7142857 1.0000000

第五步 使用UPGMA方法或者Ward法，对群落进行聚类

1

plot.hc <- hclust(plot.dist , "ave")

绘图

1
2

plot(plot.hc)
plot(plot.hc, hang=-1)

转换为dendrogram之后绘图

1
2

dendro.plot<-as.dendrogram(plot.hc)
plot(dendro.plot, horiz=T)

生物多样性指数的计算在R中计算十分简单，这里介绍一下用vegan包计算生物多样性指数的方法。

首先要将R软件和用到的程序包安装好。在数据处理和格式转换过程中，需要用到openxlsx、reshape2和vegan包。如果没有安装这几个程序包，就在R console中输入以下命令安装。

1
2
3

install.packages('reshape2')
install.packages('openxlsx')
install.packages('vegan')

第一步 在Excel中输入数据

建议将数据整理成样方、物种、个体数格式，这也是样方调查时最常见的一种格式，将数据输入到Excel中。注意:Excel文件列的名字要与下面的格式完全一致，因为R对大小写敏感，所以大小写也要一致。

1
2
3
4
5
6
7
8
9
10
11
12
13

plotname species abundance
plot1     sp1         3
plot1     sp2         6
plot1     sp3         1
plot1     sp4         2
plot1     sp5         1
plot2     sp1         8
plot2     sp3        30
plot3     sp4         2
plot3     sp2         1
plot3     sp6         1
plot3     sp7         3
.....

保存为 herbplots.xlsx，假设保存在“D:/herb/”。

第二步 读取Excel文件

1
2
3
4

library(openxlsx)
library(reshape2)
library(vegan)
herb.data <- read.xlsx("D:/herb/herbplots.xlsx")

第三步 转换为样方-物种矩阵

1
2
3
4
5

library(reshape2)
herb.mat <- acast(herb.data, 
             formula = plotname ~ species , 
             value.var = "abundance", 
             fill = 0)

此时生成的矩阵，格式如下：

1
2
3
4

       sp1 sp2 sp3 sp4 sp5 sp6 sp7
plot1   3   6   1   2   1   0   0
plot2   8   0  30   0   0   0   0
plot3   0   1   0   2   0   1   3

第四步 计算生物多样性指数

Shannon-Wiener指数

1

Shannon.Wiener <- diversity(herb.mat, index = "shannon")

Simpson指数

1

Simpson <- diversity(herb.mat, index = "simpson")

Inverse Simpson指数

1

Inverse.Simpson <- diversity(herb.mat, index = "inv")

物种累计数

1
2

S <- specnumber(herb.mat)
plot(S)

#Pielou均匀度指数

1

J <- Shannon.Wiener/log(S)

金缕梅原产于我国亚热带地区。在湖北神农架山区，干旱的山顶上，生长有大量的金缕梅。北京的金缕梅均为栽培。

金缕梅是春季最早开花的植物之一。金缕梅先花后叶，花瓣金色，狭长，在风中微微飘动，十分惹人喜爱。

现将今天拍摄的几张照片粘贴上来，与读者共享。

BIOMOD是R软件的一个程序包，是BIOdiversity MODelling的缩写，主要用来模型物种分布区及其随气候的变化。该软件是Wilfried Thuiller等三位学者在法国蒙彼利埃CNRS功能与进化生态学中心开发的。W. Thuiller是活跃在学术一线的年轻科学家，其发表的学术论文具有较大影响。

BIOMOD程序包可以运行以下模型，对物种潜在分布区进行估算，同时也能够对相应的气候变化的情形下，物种分布区的变化给出预测。预测过程中可以设置物种迁移速率，从而更接近真实结果。

• Generalised Linear Models (GLM)
• Generalised Additive Models (GAM)
• Classication Tree Analysis (CTA)
• Articial Neural Networks (ANN)
• Surface Range Envelope (SRE)
• Generalised Boosting Model (GBM)
• Breiman and Cutler’s random forest for classication and regression (randomForest)
• Mixture Discriminant Analysis (MDA)
• Multiple Adaptive Regression Splines (MARS)

希望本文档对刚刚开始使用BIOMOD的程序包的人员有所帮助。

由于译者水平有限，在翻译过程中难免有不妥的地方，还请读者及时批评指正, 若有疑问可发邮件至jinlongzhang01@gmail.com 。

译 者

2010年1月

于北京 香山