R语言计算β多样性指数及分析

Posted 2023-05-18

参考技术A 计算β多样性指数需要用到phyloseq包。它的安装方式不同于简单的install.packages（“phyloseq”）

有两种方法可以安装

1.先安装BiocManager

install.packages("BiocManager")

library("BiocManager")

BiocManager::install("phyloseq")

library("phyloseq")

2.source("https://bioconductor.org/biocLite.R")

biocLite("phyloseq")

#安装phyloseq

library("phyloseq")

安装并加载了phyloseq包后，开始读取数据，前面计算α多样性，用到的是read.table……

qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")

#读取数据，参数都是文件名，注意加后缀

#otufilename指定out表格，mapfilename指定map文件（分组数据)

#treefilename指定有根进化树文件

#refseqfilename指定代表序列文件

otu<-qiimedata@otu_table@.Data

#从qiimedata里面提取otu

sum_of_otus<-colSums(t(otu))

#t_转置,colsums计算列的和,即计算各个otu检测到的总序列数，为了筛掉一些总序列数过低的otu（可能是测序错误）

sum_of_otus

#查看otu总序列数

selected_otu<-names(sum_of_otus)[sum_of_otus>10]

#获取总序列数大于10的otu id

sub_qiimedata <- prune_taxa(selected_otu, qiimedata)

#筛选总序列数大于10的otu的phyloseq数据

weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')

#计算样本间加权unifrac

unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')

#计算样本间非加权unifrac

bray_curtis <- distance(sub_qiimedata, method='bray')

write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)

#保存距离矩阵

#计算样本间Bray-Curtis距离矩阵，method 可选" wunifrac ", " unifrac " ，"jaccard"等

pcoa_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "PCoA")

#基于Bray-Curtis距离矩阵的PCoA排序分析

p<-plot_ordination(sub_qiimedata, pcoa_of_bray_curtis, type="samples", color="Group1",shape = "Group1")

#将PCoA排序分析结果可视化

library("ggplot2")

p<-p+ scale_colour_manual(values=c("#DC143C","#808000","#00CED1")) + geom_point(size=2) +ggtitle("PCoA of Bray-Curtis distance")+theme(text = element_text(size = 15))

#修改图形大小,ggtitle加标题,stat_ellipse加椭圆

#用scale_colour_manual(values=c())自定义颜色，可查颜色的16进制对照表

p

nmds_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "NMDS")

#基于Bray-Curtis距离矩阵的NMDS排序分析

p1<-plot_ordination(qiimedata, nmds_of_bray_curtis, type="samples", color="Group1")

#将NMDS排序分析结果可视化

# color=“Group1”指定不同分组的点染不同颜色

p1

p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))

#对图片进行适当修饰， stat_ellipse()加椭圆， ggtitle()加标题

ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width

PCoA中的两个点距离，接近β多样性指数

PCA(Principal Components Analysis)即主成分分析，也称主分量分析或主成分回归分析法，首先利用线性变换，将数据变换到一个新的坐标系统中;然后再利用降维的思想，使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上，第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数，同时还保持数据集的对方差贡献最大的特征，最终使数据直观呈现在二维坐标系。

PCoA(Principal Co-ordinates Analysis)分析即主坐标分析，可呈现研究数据相似性或差异性的可视化坐标，是一种非约束性的数据降维分析方法，可用来研究样本群落组成的相似性或相异性。它与PCA类似，通过一系列的特征值和特征向量进行排序后，选择主要排在前几位的特征值，找到距离矩阵中最主要的坐标，结果是数据矩阵的一个旋转，它没有改变样本点之间的相互位置关系，只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分，而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

NMDS图中两个点的距离的排序，接近β多样性指数的排序

R语言基础题及答案——R语言与统计分析第五章课后习题(汤银才)

R语言与统计分析第五章课后习题(汤银才)

题-1

设总体$X$是用无线电测距仪测量距离的误差, 它服从$(\alpha ,\beta )$上的均匀分布, 在200次测量中, 误差为$X_i$的次数有$n_i$次:

$X_i$	3	5	7	9	11	13	15	17	19	21
$n_i$	21	16	15	26	22	14	21	22	18	25

求$(\alpha ,\beta )$的矩法估计值(注: 这里的测量误差为$X_i$是指测量误差在$(X_i-1,Xi+1)$间的代表值.

# 构造X_i,n_i序列
X_i<-seq(3,21,by=2)
n_i<-c(21,16,15,26,22,14,21,22,18,25)

# 依据X_i,n_i序列还原X序列
X<-rep(n_i,X_i)

# 均值和标准差
mu<-mean(X)
sigma<-sd(X)

# 求\\alpha和\\beata
# (a+b)/2=E(x),(b-a)^2/12=D(x)
# a+b=2E(x),b-a=2sqrt(3)sd(x)
# a=E(x)-sqrt(3)sd(X)
# b=E(x)+sqrt(3)sd(X)
alpha<-mu-sqrt(3)*sigma; print(alpha)
beata<-mu+sqrt(3)*sigma; print(beata)

[1] 13.88
[1] 27.15

题-2

为检验某自来水消毒设备的效果，现从消毒后的水中随机抽取50L，化验每升水中大肠杆菌的个数(假设1L水中大肠杆菌个数服从泊松分布), 其化验结果如下：

大肠杆菌数/$L$	0	1	2	3	4	5	6
水的升数	17	20	10	2	1	0	0

试问平均每升水中大肠杆菌个数为多少时，才能使上述情况的概率达到最大

# NUM:大肠杆菌数, v:水的升数
NUM<-0:6;
v<-c(17,20,10,2,1,0,0)

# 求期望的均值
E<-mean(NUM*v); print(E)

[1] 7.143

题-3

已知某种木材的横纹抗压力服从$N(\mu ,{\sigma }^2)$, 现对十个试件作横纹抗压力试验,得数据如下$(kg/cm2)$

$$482,493,457,471,510,446,435,418,394,469$$

1. 求$\mu$的置信水平为0.95的置信区间
2. 求$\sigma$的置信水平为0.90的置信区间

# 压力值
P<-c(482,493,457,471,510,466,435,418,394,469)

# \\mu 95%置信区间
t.test(P)$conf.int


# \\sigma 90%置信区间
# 课本上提供的函数chisq.var.test()
chisq.var.test <- function (x,var,alpha,alternative="two.sided"){
  options(digits=4)
  result<-list( )
  n<-length(x)
  v<-var(x)
  result$var<-v
  chi2<-(n-1)*v/var
  result$chi2<-chi2
  p<-pchisq(chi2,n-1)
  if(alternative == "less"|alternative=="greater"){
    result$p.value<-p
  } else if (alternative=="two.sided") {
    if(p>.5)
    p<-1-p
    p<-2*p
    result$p.value<-p
  } else return("your input is wrong")
    result$conf.int<-c(
      (n-1)*v/qchisq(alpha/2, df=n-1, lower.tail=FALSE),
      (n-1)*v/qchisq(alpha/2, df=n-1, lower.tail=TRUE))
    result
}


chisq.var.test(P,var(P),0.1)$conf.int

# =============================================+
# 以下为测试:                                  |
# ---------------------------------------------+
# x<-c(175,176,173,175,174,173,173,176,173,179)|
# t.test(x)$conf.int                           |
# ---------------------------------------------+         
# chisq.var.test(x,var(x),0.05)$conf.int       |
# ---------------------------------------------+

[1] 434.4 484.6
attr(,“conf.level”)
[1] 0.95
[1] 653.9 3327.0

题-4

某卷烟厂生产两种卷烟A和B，现分别对两种香烟的尼古丁含量进行6次试验，结果如下：

卷烟A	25	28	23	26	29	22
卷烟B	28	23	30	35	21	27

若香烟的尼古丁含量服从正态分布，

1. 问两种卷烟中尼古丁含量的方差是否相等?

2. 试求两种香烟的尼古丁平均含量差的95%置信区间?

# 香烟A\\B数据
A<-c(25,28,23,26,29,22)
B<-c(28,23,30,35,21,27)

# var test:
# p-value>0.05方差相等
var.test(A,B)

# 平均含量差的95%
t.test(x,y,var.equal = TRUE)$conf.int

F test to compare two variances
.
data: A and B
F = 0.3, num df = 5, denom df = 5, p-value = 0.2
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.04187 2.13821
sample estimates:
ratio of variances
0.2992
.
[1] -0.0007721 0.0084388
attr(,“conf.level”)
[1] 0.95

题-5

比较两个小麦品种的产量, 选择22块条件相似地试验田, 采用相同的耕作方法做实验, 结果播种甲品种的12块实验田的单位面积产量和播种乙品种的12块实验田的单位面积产量分别为：

甲品种	628	583	510	554	612	523	530	615	573	603	334	564
乙品种	535	433	398	470	567	480	498	560	503	426	338	547

假定每个品种的单位面积产量均服从正态分布,甲品种产量的方差为2140, 乙品种产量的方差为3250, 试求这两个品种平均面积产量差的置信水平为0.95的置信上限和置信水平为0.90的置信下限.

# 甲:A, 乙:B 
X_A<-c(628,583,510,554,612,523,530,615,573,603,334,564)
X_B<-c(535,433,398,470,567,480,498,560,503,426,338,547)

sigma_A<-2140
sigma_B<-3250

# 课本提供的函数two.sample.ci()
two.sample.ci=function(x,y,conf.level=0.95,sigma1,sigma2)
{options(digits=4)
  m=length(x);n=length(y)
  xbar=mean(x)-mean(y)
  alpha=1-conf.level
  zstar=qnorm(1-alpha/2)*(sigma1/m+sigma2/n)^(1/2)
  xbar+c(-zstar,+zstar)
}

# 置信水平为0.95的置信上限和置信水平为0.90的置信下限       
two.sample.ci(X_A,X_B,conf.level=0.95,sigma_A,sigma_B)[2]
two.sample.ci(X_A,X_B,conf.level=0.90,sigma_A,sigma_B)[1]

[1] 114.4
[1] 37.97

题-6

有两台机床生产同一型号的滚珠，根据以往经验知，这两台机床生产的滚珠直径都服从正态分布. 现分别从这两台机床生产的滚珠中随机地抽取7个和9个，测得它们的直径如下(单位: mmm)

机床甲	15.2	14.5	15.5	14.8	15.1	15.6	14.7
机床乙	15.2	15.0	14.8	15.2	15.0	14.9	15.1	14.8	15.3

试问机床乙生产的滚珠的方差是否比机床甲生产的滚珠直径的方差小?

# 甲:x, 乙:y 
x=c(15.2,14.5,15.5,14.8,15.1,15.6,14.7)
y=c(15.2,15.0,14.8,15.2,15,14.9,15.1,14.8,15.3)
var.test(x,y)
# ratio of variances 5.216
# E{σx^2/σy^2}=5.216
# 乙方差小于甲

F test to compare two variances
.
data: x and y
F = 5.2, num df = 6, denom df = 8, p-value = 0.04
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
1.121 29.208
sample estimates:
ratio of variances
5.216

题-7

某公司对本公司生产的两种自行车型号A、B的销售情况进行了了解,随机选取了400人询问他们对A、B的选择, 其中有224人喜欢A, 试求顾客中喜欢A的人数比例p的置信水平为0.99的区间估计.

binom.test(224,400,conf.level=0.99)

Exact binomial test
.
data: 224 and 400
number of successes = 224, number of trials = 400, p-value = 0.02
alternative hypothesis: true probability of success is not equal to 0.5
99 percent confidence interval:
0.4944 0.6241
sample estimates:
probability of success
0.56

题-8

某公司生产了一批新产品, 产品总体服从正态分布, 现要估计这批产品的平均重量, 最大允许误差为1, 样本标准差s =10, 试问在0.95的置信度下至少要抽取多少个产品?

# 课本上size.norm2( )的定义
size.norm2<-function(s,alpha,d,m){
  t0<-qt(alpha/2,m,lower.tail=FALSE)
  n0<-(t0*s/d)^2
  t1<-qt(alpha/2,n0,lower.tail=FALSE)
  n1<-(t1*s/d)^2
  while(abs(n1-n0)>0.5){
    n0<-(qt(alpha/2,n1,lower.tail=FALSE)*s/d)^2
    n1<-(qt(alpha/2,n0,lower.tail=FALSE)*s/d)^2
  }
  n1
}

# 最后一项m是事先给定的一个很大的数(课本如是说)
size.norm2(10,0.05,2,1000)

[1] 98.44

题-9

根据以往的经验，船运大量玻璃器皿，损坏率不超过5%. 现要估计某船中玻璃器皿的损坏率，要求估计与真值间不超过1%, 且置信度为0.90, 那么要抽取多少样本验收可满足上述要求?

# 课本上size.bin( )的定义
size.bin=function(d,p,conf.level){
  alpha=1-conf.level
  ((qnorm(1-alpha/2))/d)^2*p*(1-p)
}

size.bin(0.01,0.05,0.90)

[1] 1285