花粉传播有关的几个公式

$d_{j,k}=\sqrt{(x_{j}-x_{k})^2 +(y_{j}-y_{k})^2}$

$f(x, y; \delta, b)=\frac{ b\Gamma(3/b)}{2\pi\delta^2\Gamma(2/b)^3}e^{(\frac{\Gamma(3/b)d}{\Gamma(2/b)\delta})^2}$

$\pi_{r, jk}=\frac{F_{r,k}f(d_{jk}; \delta, b))}{\sum_{l:fathers}F_{r,l}f(d_{jl};\delta,b)}$

$L(g|\alpha_{COV}, \delta, b, s, m) = \prod_{o:offspring}[sT(g_{0}|g_{j_{0}}, g_{j_{0}}) + mT(g_{0}|g_{j_{0}}, BAF) +(1-s-m)\sum_{k:father}\pi_{j_{0}k}T(g_{0}|g_{j_{0}}, g_{k})]$

Putty退出后保持背景程序运行

Putty远程登录Linux，在putty会话窗口关闭后， 远程的命令也停止了执行， 如果想要命令继续执行， 则需要用screen程序。

方法如下：

  1. 用Putty登录服务器;用户名， 密码
  2. 输入screen, 给出帮助提示， 按空格跳过
  3. 此时， 程序又回到命令提示符， 但是操作都将在screen中记录。 运行程序（可能需要长时间）。 按住Ctrl + a + d， 此时， 窗口程序与Putty的会话断开， 并提示： 
  4. 1
    2
    3
    There is a screen on:
    27469.ttt-uuu (06/14/2012 10:10:58 AM) (Detached)
    1 Socket in /var/run/screen/S-helixcn.

此时关闭putty的窗口对运算结果将没有影响。

  1. 再次用putty远程登录linux时, 运行 screen -ls 显示所有正在运行的窗口。
  2. 欲回到任何之前的对话窗口， 只需输入 screen -r 27469 即可 （这里27469在每次会话中都是变化的）， 进入窗口后， 如果要结束会话， 从相应程序退出后， 输入exit退出该窗口， 从而返回到最开始登录的状态。

R读取Excel文件

R中现在一般用 openxlsx 程序包的 read.xlsx 或 write.xlsx 读写xlsx文件。

也可以使用readxl程序包的 read_excel() 函数读取 xlsx, xls文件。readxl 只能读取Excel文件，不能保存。

2018年2月20日更新

群落密度制约：寻找圆盘和圆环内的点

在基于个体水平的样地数据分析中，有时候需要找到某一个个体周围的所有个体。特别是在进行密度制约相应研究时，寻找这样的点就更为重要。

为了减少运算量，我们需要先找出某一个点上下左右（各变化所给半径后）所围成的正方形内的所有物种，再基于此小数据集，计算每个个体周围的个体数量。下面是给出两个函数，分别计算从每个个体开始，所给半径的圆内，所有出现个体的ID，以及所给圆环内，所有出现个体的ID.

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
## Individuals within the round plate given the radius from each individual

round.id <- function (dat, r = 20){
    gx <- dat$gx
    gy <- dat$gy
    id.list <- list()    for (i in 1:nrow(dat)){
        x <- gx[i]
        y <- gy[i]
        dat.sub <- dat[(gx < x + r) & 
                         (gx > x - r) & 
                         (gy < y + r) & 
                         (gy > y - r), ]        
          if(nrow(dat.sub) > 1){
            sp <- c()            
            for(j in 1:nrow(dat.sub)){
                mat <- data.frame(c(x, dat.sub$gx[j]), 
                                  c(y, dat.sub$gy[j]))
                sp[j] <- (dist(mat) <= r)
            }
            id.include <- dat.sub$ID[sp]
        } else {
            id.include <- 0
        }
        id.list[[i]] <- id.include
    }    
    return(id.list)
}

#### 

Ringring.id <- function (dat, r = c(20, 25)){
    r1 <- max(r)
    r2 <- min(r)
    gx <- dat$gx
    gy <- dat$gy
    id.list <- list()    
    
    for (i in 1:nrow(dat)){
        x <- gx[i]
        y <- gy[i]
        dat.sub <- dat[(gx < x + r1) &
                         (gx > x - r1) &
                         (gy < y + r1) &
                         (gy > y - r1), ]        
                         
           if(nrow(dat.sub) > 1){
           sp <- c()           
           for(j in 1:nrow(dat.sub)){
               mat <- data.frame(c(x, dat.sub$gx[j]), 
                                 c(y, dat.sub$gy[j]))
               sp[j] <- (dist(mat) <= r1 & dist(mat) >= r2)
           }
           id.include <- dat.sub$ID[sp]
        } else {
           id.include <- 0
        }
        id.list[[i]] <- id.include
    }    
    return(id.list)
}
### Example
ID <- 1:300
gx <- runif(300, 0, 100)
gy <- runif(300, 0, 100)
dat <- data.frame(ID, gx, gy)
res.round <- round.id(dat, r = 25)
par(mfrow = c(1, 2)) 
### Circle
plot(gy ~ gx, data = dat, col = "gray" )
vt = seq(-pi,pi,by=0.002);
x = 25 * sin(vt) + gx[10];
y = 25 * cos(vt) + gy[10];

lines(y ~ x);

points(gx[10] , gy[10] , pch = 19, col = 2)
round10 <- res.round[[10]]
points(gx[round10], gy[round10], col = "green")

##### Ring

res.ring <- ring.id(dat, c(15, 25))
plot(gy ~ gx, data = dat, col = "gray" )
vt = seq(-pi,pi,by=0.002);
x1 = 25 * sin(vt) + gx[10];
y1 = 25 * cos(vt) + gy[10];

lines(y1 ~ x1);
x2 = 15 * sin(vt) + gx[10];
y2 = 15 * cos(vt) + gy[10];

lines(y2 ~ x2);
points(gx[10] , gy[10] , pch = 19, col = 2)
ring10 <- res.ring[[10]]

points(gx[ring10], gy[ring10], col = "blue")

北欧名城乌普萨拉

北欧名城乌普萨拉 Uppsala ，因其拥有的乌普萨拉大学而闻名于世。乌普萨拉大学在世界大学的排名中名列达到83位。乌普萨拉中心城区完好得保留着古建筑，风格高贵，典雅。历史、艺术、科学，就在这座小城交织，演奏出动人的旋律。

几百年来，乌普萨拉诞生了解剖学之父鲁德贝格Olaus Rudbeckius，分类学之父林奈Carolus Linnaeus ，天文学家摄尔休斯Anders Celsius（摄氏温度的提出者），现代化学命名体系的提出者贝采利乌斯Jöns Jakob Berzelius，光谱学的奠基人埃斯特罗姆 Anders Jonas Ångström（光谱学单位Å埃为纪念他）， 物理化学大师阿伦尼乌斯Svante Arrhenius等众多的科学家。这里还养育了大量其它领域的学者以及艺术家。不愧为北欧的一颗明珠。

由于在乌普萨拉停留短暂，仅能拍一些照片，与大家分享。

于默奥的蓝天白云

于默奥位于瑞典北部，这里的空气洁净，天空湛蓝，夏天阳光明媚。于默奥是一座大学城，有瑞典农业大学和于默奥大学，近些年，这两所学校的学术研究取得了长足进步。

MCMC中的Metropolis Hastings抽样法

Metropolis Hasting（下面简称MH）是蒙特卡罗马尔科夫链中一种重要的抽样方法。本文简要介绍MH算法，并给出一个实例。

MH算法在参数空间随机取值，作为起始点。按照参数的概率分布生成随机的参数，按照这一系列参数的组合，计算当前点的概率密度。依据当前点和起始点概率密度比值是否大于(0,1)之间的随机数来判断是否保留当前点。若当前点的概率密度大于该随机数，就称这种状态为接受状态，此时，在满足参数概率分布的前提下，继续随机抽取参数的组合，作为下一点，计算下一点的概率密度，并计算下一点概率密度和概率密度的比值，并继续循环。若当前点不能被接受，则继续在满足参数概率分布的前提下，继续生成随机数，作为新的参数组合，直到参数组合能够被接受为止。

下面是用MCMC MH抽样法生成满足一定条件二元正态分布的例子，供感兴趣的同仁参考。
问题：对于一个二元正态分布，给定输入数据点向量x （x包含两个元素，x1，x2）， 平均值参数向量mu（x1，x2），和sigma（方差矩阵），写出二元正态分布的概率密度函数。

该问题引自： http://users.aims.ac.za/~ioana/

由于本人对MCMC理解不深，对一些概念的理解难免出现错误，望读者能够批评得阅读，若发现问题，请及时告知本人。

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
## 解答：
mvdnorm <- function(x, mu, sigma){
    #从x减去mu
    x.minus.mu <- x - mu 
    exp.arg    <- -0.5 * sum(x.minus.mu * solve(sigma, x.minus.mu))  
    # det(sigma) sigma 的行列式
    return( 1 / (2 * pi * sqrt(det(sigma))) * exp(exp.arg) )
}

## 问题二
## 假设二元正态分布的参数如下： 
## 两个维度的平均值分别为 2， 3
# 协方差矩阵为
# 4 1
# 1 4
# 尝试用蒙特卡洛马尔科夫链 Metropolis Hastings 抽样法生成后验分布，进行10000次随机抽样，并计算随机点的接受率。
# 答： 按照题意，有
mu <- c(2 ,3)
sigma <- matrix(c(4, 1, 1, 4), nrow = 2)
# 限制sampler在空间的移动速率，数值越大，变化越快，该数值的设定待进一步讨论。
sd.proposal <- 2
## 设定模拟的次数
n <- 10000
## 生成NA组成的矩阵，用于保存模拟的结果
x <- matrix(nrow = n, ncol = 2)
# 设定sampler的初始值，假定数据点从 0, 0开始 （实际上该sampler可以从任意点开始移动）
cur.x <- c(0, 0)
# 计算给定初始值时的概率密度
cur.f <- mvdnorm(cur.x, mu, sigma)
### 蒙特卡洛马尔科夫链
n.accepted <- 0
for(i in 1:n){
    new.x <- cur.x + sd.proposal * rnorm(2)  ## 随机生成x
    new.f <- mvdnorm(new.x, mu, sigma)       ## 计算概率密度
    if(runif(1) < new.f/cur.f){  
        ## new.f/cur.f 概率密度的比率 和 (0,1)之间的随机数相比
        ## 若该比率小于随机数，则接受该点
       n.accepted <- n.accepted + 1
       cur.x <- new.x
       cur.f <- new.f
       }
    x[i,] <- cur.x                            ## 将cur.x存到第i行
}
#查看接受率
n.accepted/n
#查看每个变量的随机变化情况
par(mfrow=c(2,1))
plot(x[,1], type="l", xlab="t", ylab="X_1", main="Sample path of X_1")
plot(x[,2], type="l", xlab="t", ylab="X_2", main="Sample path of X_2")
## 绘制椭圆概率密度图
library(MASS)
proline.density <- kde2d(x[,1], x[,2], h = 5)
par(mfrow = c(1, 1))
plot(x, col = "gray")
contour(proline.density, add = TRUE)
points(2,3, pch = 3)

抽样后的概率密度
Trace plot

将比对好的fasta序列转换成relaxed phylip格式

在推断进化树的过程中,我们有时候希望phylip格式的文件能够将序列名称完整保留下来。这里就需要用relaxed phylip格式。

不过，通常情况下，如果序列名超过10字符，ClustalX和MUSCLE等软件给出的phylip文件会将名称截取到十个字符，也就是遵循传统的PHYLIP格式。

为了将比对好的FASTA文件直接转换为 relaxed phylip格式的文件, 这里提供了R的源代码。供有兴趣的同仁参考。

请参考最新版本的 phylotools 程序包

1
2
3
4
5
library(devtools)
install_github("helixcn/phylotools")
library(phylotools)
dat <- read.fasta("aligned_fasta.fasta")
dat2phylip(dat, outfile = "out.phy")

什么是bootstrap？

Bootstrap又称自展法，是用小样本估计总体值的一种非参数方法，在进化和生态学研究中应用十分广泛。例如进化树分化节点的自展支持率等。

Bootstrap的思想，是生成一系列bootstrap伪样本，每个样本是初始数据有放回抽样。通过对伪样本的计算，获得统计量的分布。例如，要进行1000次bootstrap，求平均值的置信区间，可以对每个伪样本计算平均值。这样就获得了1000个平均值。对着1000个平均值的分位数进行计算， 即可获得置信区间。已经证明，在初始样本足够大的情况下，bootstrap抽样能够无偏得接近总体的分布。

下面是一个实例：

例如，假设有一批产品，随机抽出30个，使用寿命（天数）如下，试用bootstrap的方法估计这批产品寿命95%的置信区间。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
dat <- c(119,120,131,209,210,337,332,287,146,1
       29,232,169,208,253,142,105,419,179,
         324,287,115,132,308,356,286,221,204,
         105,45,245)

### 查看原始数据的频数直方图
hist(dat, col = "gray")
#生成一个存储器
boot.sample <- list()
## 循环1000次，有放回的抽样，每次生成的
## 新样本存储在boot.sample中
for(i in 1:1000){
     boot.sample[[i]] <- sample(dat,size = 30, replace = TRUE)
}
## 求每个样本的mean,结果为1000个bootstrap样本的mean
boot.mean <- unlist(lapply(boot.sample, mean))
## 频数直方图
hist(boot.mean, col = "gray")
## 求95%的置信区间
CI95 <- quantile(boot.mean, probs = c(0.025, 0.975))
## 在频数直方图上加置信区间
abline(v = CI95, col = "red")