什么是系统发育独立差（Phylogenetic Independent Contrast）？

系统发育独立差（Phylogenetic Independent Contrast, PIC）是去除性状分析中物种系统发育关系的一种方法，是美国进化生物学家 J. Felsenstein于1985年提出的。

当时，有动物学家检验了若干哺乳动物脑容量和体重的相关性。然而，所用统计方法都假设各个种的性状是相互独立的，采集自同一个正态分布的总体。由于不同种之间存在系统发育上的联系，因此不能直接套用常规统计模型，所以所得结果也就不可信了。在Felsenstein之前，一些学者已经意识到这个问题，但是都没能给出理想的解决方案。

Felsenstain认真思考了这个问题，并且认为：研究多个不同性状的相关性时，必须考虑系统发育关系。他提出：假设物种性状的进化服从布朗运动，那么用系统发育关系上最近的分类单元性状的差值，再经过枝长加权，就可以得到一组新的数据。这组新的数据去掉了进化信息的干扰，因此可以用常规统计方法检验。上述方法新获得的数据就称为系统发育独立差 （The Phylogenetic Independent Contrast, PIC）。

img

图1. Felsenstein (1985) 的论文

举例来说，如果要研究哺乳动物脑容量和体重的关系，就要先获得这些种之间的进化关系，也就是经过分子钟校订的进化树。然后再计算脑容量的PIC和体重的PIC，再用脑容量的PIC和体重的PIC建立线性模型，进行相关性检验。

如果假设每个性状的进化服从布朗运动模型，也就是每个子代的性状的值是从祖先节点的性状随机变化而来的，那么进化树上相邻节点的两个种性状偏离于最近共同祖先性状的期望应该是0，且二者之间的差异正比于分化时间，也就是进化树的枝长。

在1985年的论文中，Felsenstein提出了如下算法：

  1. 假设相邻的分类单元为i和j，其性状分别为X_i和X_j，i和j的共同祖先是节点为k。
  2. 先计算X_i - X_j，作为X和Y的contrast。该contrast的期望为0，方差与枝长v_i+v_j成正比
  3. 从进化树中去掉节点i和j，保留节点k，将k当做末端节点。具体来说，就是以X_i和X_j的加权平均值作为X_k，公式如下：

img

  1. 将k节点下的枝长从v_k延长为v_k + v_i**v_j*/(v_i + v_j)。这样做的原因是，第3步中，k节点的X_k并非真实值，而是用X_i 、 X_j计算到的。

经过1-4步，进化树上的末端分类单元就减少了1个，同时，也计算出了一个contrast的值。

重复1-4步，直到进化树中只有一个节点。如果进化树中有n个末端节点，那么上述过程将得到*n-*1个独立差(contrast)。图2-3展示了只有五个物种的进化树，PIC的计算公式。

img

图 2. 拥有5个分类单元的进化树，分类单元名称(1-5为末端节点，6-8为内部节点)、枝长（V）

img

图3. PIC的计算公式

ape程序包提供了计算PIC的函数pic()，该函数使用的正是Felsenstein （1985）的算法。一般用户并不需要了解函数内部的机制，只需要提供进化树和性状向量即可。

以下为ape程序包中pic函数的示例：

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
library(ape)
### The example in Phylip 3.5c (originally from Lynch 1991)
cat("((((Homo:0.21,Pongo:0.21):0.28,",   
"Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);",   
file = "ex.tre", sep = "\n")
tree.primates <- read.tree("ex.tre")
X <- c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968)
Y <- c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259)
names(X) <- names(Y) <- c("Homo", "Pongo", "Macaca", "Ateles", "Galago")
pic.X <- pic(X, tree.primates)
pic.Y <- pic(Y, tree.primates)
cor.test(X, Y)
####  Pearson's product-moment correlation
#### data:  X and Y
## t = 2.5736, df = 3, p-value = 0.08224
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1964310  0.9884173
## sample estimates:
##       cor
## 0.8296107
cor.test(pic.X, pic.Y)
####  Pearson's product-moment correlation
#### data:  pic.X and pic.Y
## t = -0.85623, df = 2, p-value = 0.4821
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9874751  0.8823934
## sample estimates:
##        cor
## -0.5179156
lm(pic.Y ~ pic.X - 1) # both regressions
#### Call:
## lm(formula = pic.Y ~ pic.X - 1)
#### Coefficients:
##  pic.X
## 0.4319
lm(pic.X ~ pic.Y - 1) 
# through the origin
#### Call:
## lm(formula = pic.X ~ pic.Y - 1)
#### Coefficients:
## pic.Y
## 0.998

从以上结果可以看出，将性状直接进行相关分析以及用性状的PIC进行相关分析，所得结果是有差异的：去除了系统发育关系之后，性状X和性状Y之间的相关性降低了。去除系统发育关系之前 （Pearson’s r = 0.8296, t = 2.5736, df = 3, p-value = 0.08224），去除系统发育关系之后（Pearson’s r = -0.5179156, t = -0.85623, df = 2, p-value = 0.4821 ）。

这里比较难以理解的是，计算出来的PIC总是比真实性状数据少一个，但是不同性状的PIC都是比真实数据少一个，用不同性状的PIC向量建立线性模型，所得结果，就是我们要计算的性状之间的相关性。

ape的pic函数要求不同种性状出现的顺序与在进化树中出现的顺序完全相同。为此，picante程序包提供了match.phylo.data()函数，只要用户为性状数据（无论是data.frame还是vector）提供了对应的物种名（如果是多个性状，格式应该为dataframe，其行名应该是物种名），该函数就能将性状数据按照进化树中出现的顺序排好，非常方便。

以下代码展示的就是怎样通过植物物种名（拉丁名）获得进化树，再计算植株高度的PIC，源数据来自S. Kembel博士的培训资料。

Felsenstein在1985年的论文中说，他提出的算法可计算多分枝结构的PIC，但是ape的pic函数不能处理含有多分枝结构的进化树，所以下面的代码调用了multi2di()函数，将进化树转换为二叉树。转换的方法是为每一个多分支节点随机添加一个很小的枝长，由于枝长很小，对结果的影响很小，几乎可以忽略不计。

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
26
27
28
29
30
31
32
33
34
35
setwd("...")
library(plantlist)
## This is plantlist 0.6.1.
# devtools::install_github("helixcn/plantlist")
library(V.PhyloMaker)
# devtools::install_github("jinyizju/V.PhyloMaker")
library(picante)
library(ape)
dat <- read.csv("trait_herb_114.csv", header = TRUE)
head(dat)
##                 species height
## 1     Stipa_albasiensis   48.5
## 2 Heteropappus_altaicus   15.6
## 3  Roegneria_alashanica   43.7
## 4 Saussurea_alaschanica   21.5
## 5 Astragalus_alaschanus   18.0
## 6      Panzerina_lanata   12.5
height <- dat$height
names(height) <- dat$species

# Build a phylogeny using phylo.maker
species <- gsub("_", " ", as.character(dat$species))
tab2 <- subset(TPL(gsub("_", " ", species)),              
select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
colnames(tab2) <- c("species", "genus", "family")
result1 <- phylo.maker(tab2, scenarios = c("S1"))
## [1] "Taxonomic classification not consistent between sp.list and tree."##           genus family_in_sp.list family_in_tree## 77 Hemerocallis  Xanthorrhoeaceae  Asphodelaceae
# write.tree(result1[[1]], "tree1.newick") 
# If you want to examine the tree
tree <- multi2di(result1[[1]])
# plot(tree)

# "If x has names, its values are matched to the tip labels of phy, otherwise its values are taken to be in the same order than the tip labels of phy."
data_matched <- match.phylo.data(phy = tree, data = height )
# Phylogenetic Independent Contrastheight_pic <- pic(data_matched$data, data_matched$phy)

参考资料

  • Felsenstein, J. (1985) Phylogenies and the comparative method. American Naturalist, 125, 1–15.
  • Jin, Y., & Qian, H. (2019). V. PhyloMaker: an R package that can generate very large phylogenies for vascular plants. Ecography, 42(8), 1353-1359.
  • Kembel, S. W., Cowan, P. D., Helmus, M. R., Cornwell, W. K., Morlon, H., Ackerly, D. D., … & Webb, C. O. (2010). Picante: R tools for integrating phylogenies and ecology. Bioinformatics, 26(11), 1463-1464.
  • Kembel S. 编， 张金龙 译， 在R中分析物种和系统发育多样性 http://blog.sciencenet.cn/home.php?mod=space&uid=255662&do=blog&id=1116212
  • Paradis, E., & Schliep, K. (2019). ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics, 35(3), 526-528.
  • Zhang, J.L. (2018). plantlist: Looking Up the Status of Plant Scientific Names based on The Plant List Database. R package version 0.6.1. https://github.com/helixcn/plantlist/

再忆铁峰

再忆铁峰

又到过年了，又想起了铁峰兄。铁峰姓朴，朝鲜族，在黑龙江齐齐哈尔长大，年长我一岁，离开这个世界已经五年了。

十年前，我还在读博的时候，有一次跟研究组的米湘成老师去黑龙江采集树木DNA条形码的样品，就在那次认识的铁峰。

那时候，铁峰还在韩国江原道大学读博士，但主要野外工作是在黑龙江的凉水和丰林自然保护区做的，以东北林业大学金光泽教授的实验室为基地。

到了哈尔滨，我和米老师先拜访了东林的金光泽老师，金老师、铁峰跟我们一起去采样。

我们坐着老式绿皮火车，先去凉水的实习基地。

初夏，东北的景色十分宜人。窗外，小山丘起伏十分平缓，一路上上下下。山上，草木葱葱郁郁地生长着。低地上是一片又一片农田，绿油油的，有时一眼望不到尽头。一座又一座不知道姓名的小村，镶嵌在绿色的地毯上。天空中一朵一朵的白云，懒洋洋地向前移动着。那天阳光特别好，风从车窗飕飕地吹进来，十分温暖。时间仿佛过得又快，又漫长，一个又一个小站，绥化、庆安、神树、带岭……铁峰就坐我身边，我们在火车上聊了一路，到了凉水，天色已晚。

第二天，我们要在这里九公顷的样地采集DNA条形码的样品，大家要一起钻林子了。

正是初夏，凉水的草爬子（全沟蜱 Ixodes persulcatus）很“厚”，虽然我很小心，还是被叮了三处：胸口、肚皮、腋窝。采了植物样品，回到住处，我赶紧脱下上衣，让铁峰帮看看后背上有没有草爬子，然后我帮他看。

确认身上和衣服上没有草爬子，我和铁峰一起在实习基地的楼道处理样品：把叶片装进小塑料袋，加变色硅胶，袋上写标本编号，并采集标本。这个工作有些枯燥，我和他边干活边闲聊。

那时候我没有出过国，很好奇他的留学经历，于是问：“你是朝鲜族，到了韩国，会不会觉得自己回到自己的国家？”

铁峰说：“虽说我是朝鲜族，但是一直说在中国长大，说汉语，只能说一些简单的朝鲜语。朝鲜族本身就是中国的一个民族，我自己就是中国人啊。”

我十分赞许他的话，连忙点头。

凉水样地采样结束后的第二天，我们又乘火车到伊春五营的丰林50公顷大样地采样。

丰林大洋地地势平缓，森林茂密，有很多大树。林下走起来并不困难，但是有时候一些倒木拦在前头，有的地段比较泥泞，这样的地段就有点儿困难，只能手拉手过去。

由于我们事先准备好了调查记录，只需要按照号牌找到相应的树就可以了，采样相对简单。

丰林大样地的记录里，最后要找一棵大青杨（Populus ussuriensis Kom.），但找了许久都没找到。最后，按照调查记录的坐标找到号码相同的一棵树，原来是一棵大黄柳（Salix raddeana Laksch.）。那一刻，大家都笑了起来，原来是鉴定错误。铁峰也笑得和孩子一样纯真。保护区的宋主任有点儿不好意思，也跟着一起说：“哪儿有大青杨啊？根本就是大黄柳！”

铁峰在丰林样地都工作了很久，和保护区的人员很熟，在回去的车上，保护区的工作人员开玩笑称呼他为“铁棒”，这时候他就低下头，笑一笑，也不做声。

晚上，大家一起吃饭。保护区工作人员非常热情，一直不停劝酒。我不胜酒力，被灌多了，感觉浑身发冷，见情况不妙，赶紧招呼铁峰。铁峰见我面色发白，就搀着我回住处。

路上，他跟我说：“幸好你喝多了，我正好才能逃出来，不然也肯定高了。”

夜色朦胧，我只恍惚记着摇晃的街灯，走不动路的双腿。他搀着我，需要的时候轻轻拍打我后背。

离开伊春，再见面好像就是CForBio大样地研究网络数据分析培训班，培训内容好像是植物功能性状分析。那次我给来自全国各地的几十位学员讲了用R怎样基于植物的功能性状建立聚类树，以进行功能多样性分析。中午，我拉着铁峰到香山上吃饭。他特别客气，争着要买单。

后来，听说他从韩国留学回国，去西北农林科技大学工作了。那之后我们经常在QQ上联系，聊各种问题，生活的、教学的、科研的，每次都能聊得很投机。像铁峰这样，能聊到一起，能一起说真话，能懂你的好友，十分难得。也就那期间，他在密度制约方向发了两篇很好的论文，我向他表示祝贺。

2015年过年，我用QQ发消息给他，祝他新年快乐。他的QQ发来消息说：你不知道吗？铁峰已经不在了。

我心头一惊，顿时说不出话来，赶紧打他的手机。接电话的是他姨，告诉我说：他的亲人正在西安处理他的后事……

电话里，我才知道，他是由妈妈一个人抚养长大；才知道他家人到学校后的种种遭遇，特别是工伤认定、补偿金等问题，很长时间才解决好。后来，我又和他家人通了几次电话，只能尽量安慰他的妈妈，说自己愿意提供力所能及的帮助。

铁峰离开这个世界五年了，总让人感觉这一切不是真的。一位优秀的青年生态学家，忽然之间，就走过了短暂的一生，怎能不让人悲伤。

如果有来世，我还愿意认识他，这个简单而纯粹的人。

纯粹的人，总是值得怀念。

计算进化树末端节点出现的时间及两节点共同祖先（MRCA）的年龄

目的：

  1. 建立进化树
  2. 获得进化树中每个种的延续时间
  3. 获得某两个种最近共同祖先（MRCA）距离现在的时间

代码

1
2
3
4
5
6
7
8
9
10
setwd("C:/Users/jlzhang/Desktop/age")
library(plantlist)    
# devtools::install_github("helixcn/plantlist")
## This is plantlist 0.6.1.
library(V.PhyloMaker) 
# devtools::install_github("jinyizju/V.PhyloMaker")
library(phytools)
library(dispRity)
library(openxlsx)
library(phytools)

用中文名建进化树

查询学名

注意xlsx中保存的是中文名

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
species <- read.xlsx("checklist.xlsx")
taxa.table(TPL(CTPL(species[,1])$SPECIES), "taxa_table.txt")
## Warning in CTPL(species[, 1]): Taxa: 红枝柴
##  could not be recognized.
## Note: too many rows to show, only the first few rows were printedTAXA_NAME :
## [1] "木荷"     "厚皮香"   "狗骨柴"   "杜鹃"     "薄叶山矾" "赤楠"    
## 
## SPECIES :
## [1] "Schima superba"           "Ternstroemia gymnanthera"
## [3] "Diplospora dubia"         "Rhododendron simsii"     
## [5] "Symplocos anomala"        "Syzygium buxifolium"     
## 
## SPECIES_FULL :
## [1] "Schima superba Gardn. et Champ."                 
## [2] "Ternstroemia gymnanthera (Wight et Arn.) Beddome"
## [3] "Diplospora dubia (Lindl.) Masam."                
## [4] "Rhododendron simsii Planch."                     
## [5] "Symplocos anomala Brand"                         
## [6] "Syzygium buxifolium Hook. et Arn."               
## 
## GENUS :
## [1] "Schima"       "Ternstroemia" "Diplospora"   "Rhododendron" "Symplocos"   
## [6] "Syzygium"    
## 
## GENUS_AUTHOR :
## [1] "Reinw. ex Blume" "Mutis ex L. f."  "DC."             "L."             
## [5] "Jacq."           "Gaertn."        
## 
## GENUS_CN :
## [1] "木荷属"   "厚皮香属" "狗骨柴属" "杜鹃花属" "山矾属"   "蒲桃属"  
## 
## FAMILY :
## [1] "Theaceae"         "Pentaphylacaceae" "Rubiaceae"        "Ericaceae"       
## [5] "Symplocaceae"     "Myrtaceae"       
## 
## FAMILY_CN :
## [1] "山茶科"   "五列木科" "茜草科"   "杜鹃花科" "山矾科"   "桃金娘科"
## 
## FAMILY_NUMBER :
## [1] "APGIII_334" "APGIII_330" "APGIII_350" "APGIII_344" "APGIII_335"
## [6] "APGIII_222"
## 
## ORDER :
## [1] "Ericales"    "Ericales"    "Gentianales" "Ericales"    "Ericales"   
## [6] "Myrtales"   
## 
## GROUP :
## [1] "Angiosperms" "Angiosperms" "Angiosperms" "Angiosperms" "Angiosperms"
## [6] "Angiosperms"
## 
## IUCN_CHINA :
## [1] "无危（LC）" "无危（LC）" "无危（LC）" "无危（LC）" "无危（LC）"
## [6] "无危（LC）"
## 
## ENDEMIC_TO_CHINA :
## [1] "否" "否" "否" "否" "否" "否"
## 
## PROVINTIAL_DISTRIBUTION :
## [1] "浙、赣、湘、黔、闽、台、粤、桂、琼、港"                
## [2] "皖、浙、赣、湘、鄂、川、黔、闽、粤、桂、滇、琼、港"    
## [3] "皖、苏、浙、赣、湘、川、闽、台、粤、桂、滇、琼、港"    
## [4] "皖、苏、浙、赣、湘、鄂、川、黔、闽、台、粤、桂、滇"    
## [5] "皖、苏、浙、湘、鄂、川、黔、闽、台、粤、桂、滇、藏、琼"
## [6] "皖、浙、赣、湘、黔、闽、台、粤、桂"                    
## 
## ALTITUDE :
## [1] "1000"                     "200-1400(2000-2800云南）"
## [3] "40-1500"                  "500-1200（-2500）"       
## [5] "400-3000"                 ""
##  [1] "Theaceae/Schima/Schima_superba"                        
##  [2] "Pentaphylacaceae/Ternstroemia/Ternstroemia_gymnanthera"
##  [3] "Pentaphylacaceae/Cleyera/Cleyera_japonica"             
##  [4] "Rubiaceae/Gardenia/Gardenia_jasminoides"               
##  [5] "Rubiaceae/Diplospora/Diplospora_dubia"                 
##  [6] "Ericaceae/Rhododendron/Rhododendron_simsii"            
##  [7] "Symplocaceae/Symplocos/Symplocos_anomala"              
##  [8] "Symplocaceae/Symplocos/Symplocos_sumuntia"             
##  [9] "Symplocaceae/Symplocos/Symplocos_lucida"               
## [10] "Symplocaceae/Symplocos/Symplocos_stellaris"            
## [11] "Myrtaceae/Syzygium/Syzygium_buxifolium"                
## [12] "Aquifoliaceae/Ilex/Ilex_ficoidea"                      
## [13] "Aquifoliaceae/Ilex/Ilex_rotunda"                       
## [14] "Aquifoliaceae/Ilex/Ilex_micrococca"                    
## [15] "Rutaceae/Evodia/Evodia_glabrifolia"                    
## [16] "Daphniphyllaceae/Daphniphyllum/Daphniphyllum_oldhami"  
## [17] "Lauraceae/Litsea/Litsea_elongata"                      
## [18] "Lauraceae/Litsea/Litsea_cubeba"                        
## [19] "Rosaceae/Photinia/Photinia_beauverdiana"               
## [20] "Moraceae/Ficus/Ficus_erecta"                           
## [21] "Myricaceae/Myrica/Myrica_rubra"                        
## [22] "Styracaceae/Styrax/Styrax_suberifolius"                
## [23] "Fagaceae/Castanopsis/Castanopsis_carlesii"             
## [24] "Fagaceae/Cyclobalanopsis/Cyclobalanopsis_gracilis"     
## [25] "Ebenaceae/Diospyros/Diospyros_morrisiana"              
## [26] "Elaeocarpaceae/Elaeocarpus/Elaeocarpus_japonicus"      
## [27] "Betulaceae/Carpinus/Carpinus_viminea"                  
## [28] "Pittosporaceae/Pittosporum/Pittosporum_illicioides"

读取查询到的学名

1
2
3
tab <- read.table("taxa_table.txt", sep = "/")
colnames(tab) <- c("family", "genus", "species")
tab2 <- tab[, c("species", "genus", "family")]

用V.PhyloMaker建树，速度较慢

1
2
3
4
result1 <- phylo.maker(tab2, scenarios = c("S1"))
write.tree(result1[[1]], "tree1.newick") 
# 保存到硬盘，用于用其他软件（如Figtree）查看
tree <- result1[[1]]

获得每个种延续的时间

注意，这个与进化树收录的分类单元密切相关，一般不宜用于研究。

脚本来源 http://blog.phytools.org/2013/10/finding-edge-lengths-of-all-terminal.html

1
2
3
4
n <- length(tree$tip.label)
ee <- setNames(tree$edge.length[sapply(1:n, function(x, y) which(y == x), y = tree$edge[, 2])], tree$tip.label)
ee <- as.data.frame(ee)
write.xlsx(ee, "terminal_branch_length.xlsx", row.names = TRUE)

查询两个分类单元最近共同祖先的节点名

1
2
3
tree$node.label <- paste(1:tree$Nnode + Ntip(tree))
name_node <- getMRCA(tree, tip = c("Cleyera_japonica", "Rhododendron_simsii"))
name_node_ind <- tree$node.label[name_node - Ntip(tree)]

绘制进化树

1
2
3
4
plot(tree)
axisPhylo(1, las = 1)
## see  
nodelabels(text=tree$node.label,  node=1:tree$Nnode+Ntip(tree))

img

1
2
3
plot(tree)
axisPhylo(1, las = 1)
nodelabels(paste("MRCA:", name_node), name_node, frame = "r", bg = "yellow", adj = 0)

img

查询两个种的最近共同祖先MRCA的时间

1
2
3
4
5
### From the dispRity Package
age_dat <- tree.age(tree)
age_dat[as.character(age_dat$elements) == name_node_ind, ]
##      ages elements
## 37 99.706       37

更多内容

用phylotools建立supermatrix

phylotools是用来建立DNA supermatrix矩阵的R程序包，是为了将比对好的基因文件（Fasta文件或Phylip文件）合并成超级矩阵（supermatrix）而编写的。

Supermatrix与Relaxed Phylip格式

supermatrix，就是将不同基因的片段分别比对后，再按照物种名合并形成的大型数据集，以便建立进化树用。

supermatrix一般为Relaxed Phylip格式。这是因为能够建立大型进化树的软件，如RAxML、Q-Tree等，都能识别这种格式。

Relaxed Phylip格式是Phylip格式的扩展，即第一行的两个数字分别声明序列数和位点数，两个数字中间有一个空格。第二行以序列名开始，序列名称和序列数据之间，有一个空格。传统的Phylip文件，序列名称只允许10个字符，如果超过10个字符，Phylip等软件在建树时会将名称截断，只取前10个字符。Relaxed Phylip格式则没有这种限制。

建立supermatrix的几个常用片段

为什么建立进化树要用多个DNA片段呢？这主要是因为DNA条形码的单个基因一般都太短，难以全面反映物种之间的关系，所以建立更为可信的进化树，一般需要更长的序列，即所谓多基因策略。

不同的序列，由于编码功能的差别，有些突变的速度快，有些速度慢。有些基因片段就较为稳定，在不同科属之间差异较小，而有些则变异很快，在同一个种不同个体之间都有很大的区别。在DNA条形码研究中，常用比较稳定的rbcLa基因确定植物物种的大致类群，一般可确定到科的水平；用matK等变异稍快的类群鉴定到属甚至到种的水平；用trnH-psbA等变异速率非常快的非编码基因确定更精确的种或种下的水平等。

上述三个DNA条形码序列，都是属于叶绿体基因组。当然，除上述基因外，还有很多别的基因，如ITS等也可用来建树。

DNA条形码序列一般是对PCR产物用传统的桑格测序法获得的，对样品和反应条件的要求较高，有一定的失败率。理想情况下， 一份样品的几个条形码序列都应该是完整的。但是由于各种各样的原因，有些种的某些序列总是不能正常扩增和测序，造成数据缺失。

supermatrix的基本内容

在建树开始之前，每个基因的序列需要单独比对，然后再拼接在一起。当然，这时候需要保证同一个种的不同序列要接在一起，但是位点的顺序不能错位。与此同时，若某一个种的某一个基因的序列缺失，则其对应的全部位点都需要用?补齐。

类似软件

用比对好的序列文件建立supermatrix，除了phylotools，还有以下软件：

  1. Biopython (https://biopython.org/wiki/Concatenate_nexus)
  2. Utensils (https://github.com/ballesterus/Utensils)
  3. Geneious (https://assets.geneious.com/manual/11.1/GeneiousManualsu53.html)
  4. catfasta2phyml (https://github.com/nylander/catfasta2phyml)
  5. SEDA (http://www.sing-group.org/seda/ , https://www.biostars.org/p/332853/)
  6. SuperMatrix: (http://coleoguy.blogspot.com/2015/04/r-concatenate-alignments-into.html)
  7. concat (https://figshare.com/articles/concatenation_R/3839538)
  8. AMAS (https://www.ncbi.nlm.nih.gov/pmc/articl::es/PMC4734057/)

其中Genious可以通过图形界面实现。

Phylotools的特点

  • phylotools全部用R写成，操作简单。
  • phylotools假设同一个种，在不同的fasta或者phylip文件中有相同的名称。
  • phylotools的supermatrix可以识别比对好的fasta文件或者phylip文件，生成的supermatrix以Relaxed Phylip格式保存。

phylotools的安装

phylotools的最新版保存在github上，而并未上传到CRAN。要安装phylotools，可使用如下命令：

1
devtools::install_github("helixcn/phylotools")

若未安装devtools，可使用 install.packages("devtools")安装

phylotools程序包函数列表

  • clean.fasta.name 用于清除fasta序列名中的特殊字符，以便建树的时候能够正常输入和显示
  • dat2fasta 将读取为data.frame的DNA序列转换为fasta文件
  • dat2phylip 将读取为data.frame的DNA序列转换为phylip格式
  • get.fasta.name 获取fasta文件中所有序列的名称
  • get.phylip.name 获取phylip文件中所有序列的名称
  • read.fasta 读取fasta文件
  • read.phylip 读取phylip文件
  • rename.fasta 重命名fasta文件
  • rm.sequence.fasta 从fasta文件中删除某一条特定的序列
  • split_dat 根据序列所在的分组，将数据分成独立的fasta文件，例如，可用于trnH-psbA序列按照目(orders)分组
  • sub.taxa.label 替换进化树中的物种名
  • supermat 用已经比对好的fasta文件或者phylip文件建立supermatrix超级矩阵以及RAxML分隔文件。

用phylip文件生成Relaxed Phylip格式的supermatrix

假设现有三个Phylip文件，分别编码不同的基因，建立者三个基因的supermatrix，并保存为Relaxed Phylip格式，同时生成RAxML Partition File，记录每个基因的起讫位点。

相应代码如下：

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
cat("6 22",
"seq_1 --TTACAAATTGACTTATTATA",
"seq_2 GATTACAAATTGACTTATTATA",
"seq_3 GATTACAAATTGACTTATTATA",
"seq_5 GATTACAAATTGACTTATTATA",
"seq_8 GATTACAAATTGACTTATTATA",
"seq_10---TACAAATTGAATTATTATA",file = "matk.phy", sep = "\n")

cat("5 15",
"seq_1  GATTACAAATTGACT",
"seq_3  GATTACAAATTGACT",
"seq_4  GATTACAAATTGACT",
"seq_5  GATTACAAATTGACT",
"seq_8  GATTACAAATTGACT",
file = "rbcla.phy", sep = "\n")

cat("5 50",
"seq_2 GTCTTATAAGAAAGAATAAGAAAG--AAATACAAA-------AAAAAAGA",
"seq_3 GTCTTATAAGAAAGAAATAGAAAAGTAAAAAAAAA-------AAAAAAAG",
"seq_5 GACATAAGACATAAAATAGAATACTCAATCAGAAACCAACCCATAAAAAC",
"seq_8 ATTCCAAAATAAAATACAAAAAGAAAAAACTAGAAAGTTTTTTTTCTTTG",
"seq_9 ATTCTTTGTTCTTTTTTTTCTTTAATCTTTAAATAAACCTTTTTTTTTTA",
file = "trn1.phy", sep = "\n")

supermat(infiles = c("matk.phy", "rbcla.phy", "trn1.phy"))

用fasta文件生成supermatrix

要使用比对好的fasta文件生成supermatrix， 则命令改为对应的fasta文件名即可：

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
cat(">seq_1",  "--TTACAAATTGACTTATTATA",
">seq_2",  "GATTACAAATTGACTTATTATA",
">seq_3",  "GATTACAAATTGACTTATTATA",
">seq_5",  "GATTACAAATTGACTTATTATA",
">seq_8",  "GATTACAAATTGACTTATTATA",
">seq_10", "---TACAAATTGAATTATTATA",
file = "matk.fasta", sep = "\n")

cat(
">seq_1", "GATTACAAATTGACT",
">seq_3", "GATTACAAATTGACT",
">seq_4", "GATTACAAATTGACT",
">seq_5", "GATTACAAATTGACT",
">seq_8", "GATTACAAATTGACT",
file = "rbcla.fasta", sep = "\n")

cat(
">seq_2", "GTCTTATAAGAAAGAATAAGAAAG--AAATACAAA-------AAAAAAGA",
">seq_3", "GTCTTATAAGAAAGAAATAGAAAAGTAAAAAAAAA-------AAAAAAAG",
">seq_5", "GACATAAGACATAAAATAGAATACTCAATCAGAAACCAACCCATAAAAAC",
">seq_8", "ATTCCAAAATAAAATACAAAAAGAAAAAACTAGAAAGTTTTTTTTCTTTG",
">seq_9", "ATTCTTTGTTCTTTTTTTTCTTTAATCTTTAAATAAACCTTTTTTTTTTA",
file = "trn1.fasta", sep = "\n")

supermat(infiles = c("matk.fasta", "rbcla.fasta", "trn1.fasta"))

参考文献

  • Kress, W. J., Erickson, D. L., Jones, F. A., Swenson, N. G., Perez, R., Sanjur, O., & Bermingham, E. (2009). Plant DNA barcodes and a community phylogeny of a tropical forest dynamics plot in Panama. Proceedings of the National Academy of Sciences, 106(44), 18621-18626.
  • 裴男才, 张金龙, 米湘成, 葛学军 (2011) 植物DNA条形码促进系统发育群落生态学发展, 生物多样性, 19(3):284–294
  • Roquet, C., Thuiller, W., & Lavergne, S. (2013). Building megaphylogenies for macroecology: taking up the challenge. Ecography, 36(1), 13-26.
  • https://cme.h-its.org/exelixis/web/software/raxml/index.html

怎样确定一个地区标本采集的完整程度

问题：

现有若干县的植物标本记录，即每个种在该县境内采集的标本数量，求这些县植物标本采集的完整程度。

解答：

标本采集的完整程度可以用随标本数量增加，物种累计数的快慢表示。

假设A县采集了5000份标本，而这5000份标本总共记录了500个种；B县采集了5000份标本，而这5000份标本包括了2000个种，很显然，A县的标本采集是比较充分的；相比之下，B县就很不充足，因为B县境内，物种数仍然随着标本数快速增加。

那么如何度量物种累计的快慢呢？中科院植物所的阳文静博士（Yang et al. 2013）在度量中国各县植物标本采集完整程度时，采用了物种累计数达到90%时的斜率。也就是计算标本数达到90%-100%时物种累计曲线的平均斜率。该思路简化下来，就是先计算标本数达到90%时对应的物种数，然后取平均值，然后将所有标本数对应的总物种数作为另外一个坐标点，两点确定一条直线，确定这条直线的斜率（原文信息不够详尽，可能本人理解有误）。

另一种思路是计算物种数达到90%时的切线斜率。由于是随机抽样，物种数从<90%到跨越90%的一瞬间，物种数会增加，随机抽样时，物种数会围绕某个值上下波动，因此这里可以计算多次，取平均值。

这里只给出第二种思路的R代码，即求切线的斜率。

图片

图1. 要读取的标本记录格式

图片

\2. 阳文静（Yang et al., 2013）中国各县植物标本采集完整程度

图片

图片

图3 物种数累计曲线：标本采集完整的县斜率小，标本采集不完整的县斜率大

R脚本

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
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
setwd("C:/Users/helixcn/Desktop/completeness")
library(openxlsx)

# 获得某个文件夹的所有xlsx文件字符串向量
all_files <- list.files()
excel_files <- all_files[grepl(".xlsx", all_files)]

# 针对每个标本数量，例如从30000个标本中，随机抽26000个标本，这里要随机抽取100次， nperm就是控制随机抽取的数量
nperm = 100

slope_90_100 <- NA # 初始化, slope_90_100 保存每个文件（也就是每个地点的slope值）

# k为excel文件的数量
for(k in 1:length(excel_files)){
    
    # 读取一个excel文件
    specimen_names <- read.xlsx(excel_files[k])[,1]
    
    # 求物种数
    nspecies <- length(unique(specimen_names))
    
    # 物种数不能少于10
    if(nspecies < 10){
        stop("Number of species should be more than 10")
    }
    
    # 标本数不能少于50
    if(length(specimen_names) < 50){
        stop("Number of specimens should be more than 50")
    }
    
    # 90%的物种数
    nspecies_90 <- nspecies * 0.9
    
    ## 为了减少计算量，这里只计算物种数达到90%以后的物种累计曲线
    ## 因此，先要计算在标本达到总标本数的多少时（求一个比例ratio），物种数不超过90%

    ratio = 1.0 # 设定ratio的初始值，然后猜100次，一直到标本数刚刚少于90%的物种数为止。每一步递减0.01，猜100次（nguess），足以覆盖所有的情况
    nguess = 100
    for(m in 1:nguess){
        if(length(unique(sample(specimen_names, ceiling(length(specimen_names) * ratio)))) > nspecies_90) {
             ratio <- ratio - 0.01
         } else {
             ratio <- ratio
         }
    }

    n_species_90_to_100 <- NA
    # 从能达到90%物种数的标本比例开始计算随标本数增加的物种累计曲线
    for (j in floor(length(specimen_names) * ratio):length(specimen_names)){
        n_species_90_to_100_i <- NA # 初始化随机抽取100次的平均物种数
        for (i in 1:nperm){
            n_species_90_to_100_i <- c(n_species_90_to_100_i, length(unique(sample(specimen_names, size = j))))
        }
        n_species_90_to_100 <- c(n_species_90_to_100, mean(na.omit(n_species_90_to_100_i)))  #
    }

    # 初始化的时候，NA要去掉
    n_species_90_to_100_clean <- na.omit(n_species_90_to_100)

    # 只取物种数大于百分之九十物种数时的情况
    n_species_90_to_100_clean <- n_species_90_to_100_clean[n_species_90_to_100_clean > nspecies_90]

    # 为了求得每增加一份标本，平均增加多少种，要将所得结果错位，然后相减
    par1 <- c(NA, n_species_90_to_100_clean)
    par2 <- c(n_species_90_to_100_clean, NA)

    # 相减之前，要先去掉NA，否则平均值也是NA
    # slope_90_100为一个向量，从NA开始，长度不断增加，注意slope_90_100在循环外面已经定义了
    # mean(na.omit(c(par2 - par1))) 就是从物种数达到90%之后，每增加一份标本，平均增加的物种数

    slope_90_100 <- c(slope_90_100, mean(na.omit(c(par2 - par1)))) # 求斜率，以总标本数为基准

    # 显示每个文件，slope的值
    print(paste("The average slope (90-100%) for file", excel_files[k] ,"is:", mean(na.omit(c(par2 - par1)))))
}
# 将所有文件的slope，保存到一个CSV中，同时提供文件名
write.csv(data.frame(site = excel_files, slope_90_100 = na.omit(slope_90_100)), "survey_completeness.csv")

致谢

  • 感谢中科院植物所刘慧圆老师一起讨论问题

参考文献

  • Yang, W., Ma, K., & Kreft, H. (2013). Geographical sampling bias in a large distributional database and its effects on species richness-environment models. Journal of Biogeography, 40(8), 1415–1426. https://doi.org/10.1111/jbi.12108

从DNA序列到系统发育多样性-在环保部环科院的讲座

讲座内容包括：

1 进化树及相关基本概念
2 DNA 序列的获取和比对
3 进化模型及其筛选
4 建立进化树的方法
5 启发式搜索和 bootstrap 可靠性评估
6 贝叶斯法建树
7 分化时间的校对：分子钟
8 通过植物学名建立进化树
9 系统发育多样性分析： ape、 vegan 和 picante

本讲座是2019年9月16日应环保部环境科学研究院的刘勇波副研究员（http://www.craes.cn/zjhky/craes_kydw/sssds/201808/t20180819_454216.shtml）的邀请而准备的，在此特别感谢刘师兄的邀请和热情款待。

特别感谢沈泽昊、张英、杨拓、高芳銮、向小果、刘勇波、俞方圆、曾思金等老师提出宝贵意见！

更多内容敬请关注 http://blog.sciencenet.cn/u/zjlcas 和 https://github.com/helixcn 。

已知经纬度计算两点间地理距离

1 问题

已知巴黎的地理坐标为东经2°20’14’’，北纬48°50’11’’，华盛顿的地理坐标为西经 77°03’56’’，北纬 38°55’17’’，求两城间的距离（单位千米）。

2 分析：

已知经纬度，求两点间的地理距离，不能直接用经纬度求平面距离，而是要计算球面距离（Great Circle Distance）。

如果对精度要求不高，可以用低精度公式，计算过程相对简单。如果对经度要求较高， 则要用高精度公式。高精度公式一般用于天文计算、大地测量等。在距离较远的城市之间， 两种方法所得结果可相差几十千米。

当然，现在全部通过计算机程序计算，都是直接使用高精度方法了。本文参考比利时天文学家Jean Meeus (1991)的《天文算法》81-86页的公式，给出相应的R代码。

2.1 将经纬度的度、分、秒转换为十进制

定义函数deg2dec完成这一转换:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
deg2dec <- function(h,m,s){   
   if(h < 0){
       m = - m
       s = -s
   }
   res = h + m/60 + s/3600
   return(res)
}

## 巴黎
L1 = deg2dec(-2,20,14)
phi1 = deg2dec(48, 50, 11)

## 华盛顿
L2 = deg2dec(77,03,56)
phi2 = deg2dec(38,55,17)

2.2 低精度公式

由于地球的扁率较小，可近似看作球体，利用球面三角公式即可求出距离，R代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
geodist <- function(L1, phi1, L2, phi2){
    Ri = 6371
    d = (Ri * acos(sin(phi1 * (pi / 180)) *
                      sin(phi2 * (pi / 180)) +
                      cos(phi1 * (pi / 180)) *
                      cos(phi2 * (pi / 180)) *
                      cos((L1 - L2) * (pi / 180))))
    res <- round(d, 3)    
    return(res)
}

geodist(L1, phi1, L2, phi2)
## [1] 6165.597

因此，如果假设地球是球体，华盛顿和巴黎之间的距离是6165.597km

2.3 高精度公式

假设地球是椭球，a是长轴，f是扁率，公式如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
hageodist <- function(L1, phi1, L2, phi2){
   a = 6378.14
   f = 1/298.257
   F = (phi1 + phi2)/2
   G = (phi1 - phi2)/2
   ramda <- (L1 - L2)/2
   S = (sin(G * pi / 180)^2)*(cos(ramda * pi / 180)^2) + 
       (cos(F * pi / 180)^2)*(sin(ramda * pi / 180)^2) 
   C = (cos(G * pi / 180)^2)*(cos(ramda * pi / 180)^2) + 
       (sin(F * pi / 180)^2)*(sin(ramda * pi / 180)^2) 
   omega = atan(sqrt(S / C))
   R = sqrt(S * C)/omega
   D = 2 * omega * a
   H1 = (3 * R - 1)/(2 * C)
   H2 = (3 * R + 1)/(2 * S)
   res = D * (1 + f * H1 * (sin(F * pi / 180)^2)*
               (cos(G * pi / 180)^2) -
               f * H2 * (cos(F * pi/180)^2)*
               (sin(G * pi / 180)^2))   
   return(round(res, 3))
}

hageodist(L1, phi1, L2, phi2)
## [1] 6181.628

因此，华盛顿和巴黎之间的距离是6181.628km。高精度和低精度公式的结果相差16km。

2.4 注意事项

Meeus这本书中，东经为负，西经为正，公式采用国际天文学联合会IAU1976规定的相关参数，其坐标系未指定。

注意，Meeus的算法中，经度与常用的WGS84坐标系正好相反。在WGS84坐标系中，东经为正，西经为负。现在通常所说的经纬度是用GPS获得的WGS84经纬度，所以最保险的方法是以下R脚本计算：

3 推荐的R程序包和脚本

3.1 sp包

1
2
3
4
5
6
7
8
9
10
11
12
13
library(sp)
longitude <- as.numeric(c(L1, L2))
latitude  <- as.numeric(c(phi1, phi2))
cities <- c("Paris", "Washington")
location <- cbind(longitude = -longitude, latitude)
row.names(location) <- cities
location <- data.frame(location)
coordinates(location) <- ~longitude+latitude
proj4string(location) <- CRS("+proj=longlat +datum=WGS84") # 明确使用WGS84坐标系
spDists(location)
##          [,1]     [,2]
## [1,]    0.000 6181.626
## [2,] 6181.626    0.000

3.2 geosphere包

1
2
3
4
5
6
xy <- rbind(paris =c(-L1, phi1), washington =c(-L2, phi2))
library(geosphere)
distm(xy)
##         [,1]    [,2]
## [1,]       0 6181632
## [2,] 6181632       0

4 在地图上绘制Great Circle Distance

1
2
3
4
5
6
7
8
9
10
11
library(maps)
library(geosphere)
map('world',col="grey", fill=FALSE, bg="white", lwd=0.05,interior = FALSE)
points(x = -L1, y = phi1, pch = 19)
points(x = -L2, y = phi2, pch = 19)
text(x = -L1, y = phi1 + 5, labels = "Paris", col = "blue")
text(x = -L2, y = phi2 + 5, labels = "Washington", col = "blue")

connection <- gcIntermediate(c(-L1, phi1), c(-L2, phi2), 
    n=100, addStartEnd = TRUE, breakAtDateLine=FALSE)             
lines(connection, col="red", lwd=2)

5 参考

冰糕、奶糕

华北，6月中旬，深蓝的天空万里无云。盛夏的太阳散发着威力，炙烤着金黄色的麦田。干热风吹过，麦穗沙沙作响。此时，麦穗早已又干又脆，轻轻一碰，麦粒就会掉在地上。

晌午，割麦子的人挥汗如雨。镰刀轻轻往回一收，麦杆就清脆地断开了。每割一下，地上就升起一丝尘土，仿佛一缕淡淡的青烟。每呼吸一下，干燥的空气就带走嗓子里的一点水分，仿佛要把肺抽干。塑料水壶里的水早已被烤得发烫，喝上几大口，也不觉得多解渴。

远远的，有一个小白点来了。借着麦浪，声音传得更远：“冰糕、奶糕、刨冰……”：一个不到二十岁的姑娘，推着一辆飞鸽牌二十八加重自行车，车后座挂着两个冰棍儿箱子。这是二十多年前，麦收时候的一天中午。

卖冰棍儿姑娘的出现，比凉爽的风更让人激动。

掀开冰棍儿箱的盖子，一股凉气让拿冰棍儿的手别提多舒服了。冰棍儿上面盖着一个白色的小棉被。掀开被子，冰棍儿、雪糕整齐码着：这边是奶油，那边是小豆，还有巧克力的，“俩色俩味”的……另一个箱子里是“刨冰”，坚韧的塑料袋里冻着糖水，这个“刨冰”也很受欢迎。

干渴至极的时候能吃到冰棍儿，让人难掩内心的激动。割麦子的人赶紧掏出两元钱，冰棍儿两毛，刨冰五毛，不仅仅是给一起下地干活的孩子，自己也要赶紧吃上两块消暑。

如果天气不热，大部分人是不舍得买冰棍儿吃的。只有在烈日下，冰棍儿的销路才好。今天生意确实不错。

姑娘扶了扶草帽。汗水从额头流下来，沾湿了头发，她赶忙擦了擦汗。口渴难耐，她赶紧打开水壶，灌了几口早上在家里罐好的井水。“的确良”衬衫也早已湿透。还好，还差几块钱，本钱就上来了，这样一来，剩下的二十几只冰棍儿就全是今天赚的，也许能挣上十几块。想到此处，姑娘心里宽慰了许多。

冰棍儿都是一早上从二十多里路以外的县城趸回来的。天蒙蒙亮就要出发，那还是早上四点多。

姑娘从小到大一直在学校读书：家里人都觉得读书是改变命运的唯一希望，所以一直全力支持。只是姑娘高考发挥失常，没有考上大学，落榜在家。

没有什么选择，只有安心务农，早早找户人家嫁了。出嫁之前要努力挣点儿钱，给自己攒点儿嫁妆。

之前，她没受过多少日晒雨淋，都一直是在教室里学习，听老师讲语、数、外、史、地、生，做各种练习题，参加各种测试，对未来充满憧憬，从未想过自己要去卖冰棍儿。

住校的几年，她一直都是吃家里带的咸菜，营养不良，身体瘦弱。自行车上挂上两个空冰棍儿箱子，她觉得很不好控制，特别是上下车，右腿都要跨过车的大梁，很不习惯。

从家里出发一个小时，终于骑到了冰棍儿厂了。天早已大亮，太阳升起来，照着稚嫩的脸。前面还有很多人在往箱子里装冰棍，趸冰棍儿要挨个儿。好不容易交上钱，装好箱，两个箱子忽然仿佛装满了两大桶水，沉重异常，车子就更难控制。

来到田间小路，她不敢再骑，生怕撞到树上或者拐进沟里，瘦小的她只能慢慢向前推着车子，边走，边吆喝，“冰糕、奶糕、刨冰……”。她知道，只有最暴晒的麦田里，冰棍儿才有人买。就这样，她在烈日下转了几个钟头，早已有点儿头晕目眩，中午也只是啃了两口凉馒头就咸菜。

太阳慢慢转到西边了，没那么热了，现在还有十几根冰棍儿没有卖出去。虽然天气不太热了，但是箱子里开始暖和起来，冰棍儿变得很软，马上就要化了。

她也赶紧回到村里，将软到快要拿不起来的贱卖处理，原来五毛钱一块的，现在一块钱给三块、四块。还好，最终只是还剩几块在箱子里融化了，箱子里，糖水上漂着几根木棍儿……辛苦了一整天，挣了几块钱，姑娘心里说不出的滋味。

天快黑了，她回到家。

父母没有说她。父亲沉默不语，在一旁抽着烟。母亲帮她把箱子擦干净。

晚饭很简单，炒土豆，煮凉汤（过水面）。她吃了几口，眼泪落下来，然后独自回房间了。一个高考落榜姑娘的一天就这样结束了。

其实，这里的姑娘，跟我高考落榜的姑姑所经历的，并没有什么区别。我两个姑姑在高中毕业后，在九十年代都卖过冰棍儿。二十多年后的今天，我还能想起姑姑冰棍儿没有卖完，回来后一家人一起吃冰棍儿的情景。

行走在烈日下，我经常想起麦田边上那个由远及近的小白点，一个瘦弱的姑娘，推着自行车，后面挂着冰棍儿箱子，不停吆喝着：“冰糕、奶糕、刨冰……”

描写一流科学家需要怎样的写作技巧？《张益唐：天才的野心》读后感

写文章无外乎要传递信息，本文要传播什么信息呢？那就是张益唐这样一位天才，与常人有何不同，他过着怎样的生活。作者文笔一流，用对比、烘托等笔法，写出了张益唐如神仙般的存在：一位孤独、谦逊、深邃的谦谦君子。

本文的核心思想正如标题所讲，要突出的是张益唐诗意一样的野心：专注于大问题，专注到不管庭前花开花落，达到宠辱不惊的程度——除了纯粹数学，其他都完全不在乎。

如果成功解决这些“经典问题”，数学家本人将名垂青史，但如果没有成功解决，就什么都不存在。这是一种非常危险的研究方式，绝大多数学者都不敢选择，但是张益唐这么做了，而且做到了。作者抓住了这一点，将张益唐的经历写得十分精彩。

在作者笔下，张不是刻意忽略名利与社会地位，而是完全脱离现实中的名利，只有无时无刻对自己感兴趣数学问题的不停思考。换句话说，是为了学术而学术，完全过着一种理想主义的生活。张教授以一种现代学术界无法接受的方式存在着，也就是“不按套路出牌”。

一般来说，无论你从事那一领域的研究，最少在几年时间，要有若干成果，也就是有足够数量的论文发表，否则是无法获得项目资助的，业界也认为你从此沉沦。如果几年不发表一篇论文，一般业界会认为你的职业生涯已经结束。

当代学术界的大多数论文，读者只有从事最相近研究领域的同行。而外界看一个学者的能力和重要性，多多少少都只看对方有多少第一作者的文章，看论文发表在什么刊物上，看谷歌学术的引用次数，甚至是ResearchGate分数等。在这种情况下，学术贡献本身是用一个数字衡量的，至于工作是否实现了某些方面从零到一的突破，很多情况下直接用上述“加减法”去衡量，当然是十分荒谬的，但这是目前学界广泛认同的评价体系，绝大部分人不得不接受这种评价方式。

张益唐当然是另类中的另类，假如真如文章所写，几十年只发表三篇论文，是绝对不可能在学术机构里面获得广泛认可的，因为从论文数量上看，根本不可能达到评教授的要求。在国内目前的学术环境下，也更不可能容忍，恐怕也早已被扫地出门。

本文的作者大量运用了对比的手法，突出了张益唐过往经历的跌宕起伏与与众不同。“文似看山不喜平”，将这种反差写得淋漓尽致，因此读者才愿意看下去。

为了营造氛围，同时将人物写得生动，作者同时提供了大量细节，并引用张益唐及其亲友的话，不断强化上述反差（衣食住行等等方面，都写得很详细，让人觉得可信）。最后，作者再通过学生和朋友的第三方描述，让人物的形象更丰满。

所以说，本文作者用很纯熟的笔法，把张益唐刻画得很生动、性格很鲜明，也更深入人心，从而让一位儒雅的大学者形象跃然纸上。

这也是需要很长时间训练才能达到的写作技巧。

进一步阅读

《张益唐：天才的野心》原文链接： https://mp.weixin.qq.com/s/7cCpCAwqjIYUsARWSQZrlw

一个群落beta多样性分析的完整R脚本

本文给出Abbas S. et al. (2019)的完整R脚本，内容包括群落稀疏化曲线、beta多样性和指示种分析、mantel检验、NMDS非参数排序、空间关系的PCNM转换和前向筛选、方差分解等各种分析以及作图等，供感兴趣的同行参考。

本R脚本全部为本人编写，因群落数据本人无权公开，这里只公布R脚本。篇幅所限，部分数据分析段落有删节。如您有任何问题或者建议，请跟本人联系（微信 jinlongzhang01）。

因水平有限，脚本中难免出现错误，欢迎各位读者指正。

R脚本各部分内容如下：

  • 1 稀疏化曲线，用于比较个体数不同样方的物种丰富度
  • 2 Mantel检验两个距离矩阵的相关性是否显著
  • 3 指示种分析、PCNM、环境因子的前向筛选、物种组成的方差分解以及显著性检验
  • 4 betadisper检验beta多样性各组的离散度是否有统计差异，Mantel Correlogram展示距离矩阵相关的精细结构
  • 5 分析不同演替阶段生活型组成的变化
  • 6 提取丰富度最高的种
  • 7 NMDS分析，用以展示不同演替阶段样方物种组成相似性的关系

下载完整PDF文档 res.pdf。PDF文件中部分代码语句没有换行，造成显示不完整，关注本人微信公众号ecoinformatics（二维码见本文最后），直接查看html版本的代码。

参考文献：

  • Abbas, S., Nichol, J. E., Zhang, J., & Fischer, G. A. (2019). The accumulation of species and recovery of species composition along a 70-year succession in a tropical secondary forest. Ecological Indicators, 106, 105524. DOI:10.1016/j.ecolind.2019.105524

img

img

img

img

img

img

img

img

img

img

img