有两种方法可以安装
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@[email protected]
#从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图中两个点的距离的排序，接近β多样性指数的排序
vector：向量 numeric：数值型向量 logical：逻辑型向量 character；字符型向量 list：列表 data.frame：数据框c：连接为向量或列表 sequence：等差序列 rep：重复
length：求长度 subset：求子集 seq，from:to， NA：缺失值 NULL：空对象 sort，order，unique，rev：排序 unlist：展平列表 attr，attributes：对象属性
mode，typeof：对象存储模式与类型 names：对象的名字属性
character：字符型向量 nchar：字符数 substr：取子串 format，formatC：把对象用格式转换为字符串 paste，strsplit：连接或拆分
charmatch，pmatch：字符串匹配 grep，sub，gsub：模式匹配与替换
complex，Re，Im，Mod，Arg，Conj：复数函数
factor：因子 codes：因子的编码 levels：因子的各水平的名字 nlevels：因子的水平个数 cut：把数值型对象分区间转换为因子
table：交叉频数表 split：按因子分组 aggregate：计算各数据子集的概括统计量 tapply：对“不规则”数组应用函数
+, -, *, /, ^, %%, %/%：四则运算 ceiling，floor，round，signif，trunc，zapsmall：舍入 max，min，pmax，pmin：最大最小值
range：最大值和最小值 sum，prod：向量元素和，积 cumsum，cumprod，cummax，cummin：累加、累乘 sort：排序 approx和approx fun：插值 diff：差分 sign：符号函数
abs，sqrt：绝对值，平方根 log, exp, log10, log2：对数与指数函数 sin，cos，tan，asin，acos，atan，atan2：三角函数
sinh，cosh，tanh，asinh，acosh，atanh：双曲函数
beta，lbeta，gamma，lgamma，digamma，trigamma，tetragamma，pentagamma，choose ，lchoose：与贝塔函数、伽玛函数、组合数有关的特殊函数
fft，mvfft，convolve：富利叶变换及卷积 polyroot：多项式求根 poly：正交多项式 spline，splinefun：样条差值
besselI，besselK，besselJ，besselY，gammaCody：Bessel函数 deriv：简单表达式的符号微分或算法微分
array：建立数组 matrix：生成矩阵 data.matrix：把数据框转换为数值型矩阵 lower.tri：矩阵的下三角部分 mat.or.vec：生成矩阵或向量 t：矩阵转置
cbind：把列合并为矩阵 rbind：把行合并为矩阵 diag：矩阵对角元素向量或生成对角矩阵 aperm：数组转置 nrow, ncol：计算数组的行数和列数 dim：对象的维向量
dimnames：对象的维名 row/colnames：行名或列名 %*%：矩阵乘法 crossprod：矩阵交叉乘积（内积） outer：数组外积 kronecker：数组的Kronecker积
apply：对数组的某些维应用函数 tapply：对“不规则”数组应用函数 sweep：计算数组的概括统计量 aggregate：计算数据子集的概括统计量 scale：矩阵标准化
matplot：对矩阵各列绘图 cor：相关阵或协差阵 Contrast：对照矩阵 row：矩阵的行下标集 col：求列下标集
solve：解线性方程组或求逆 eigen：矩阵的特征值分解 svd：矩阵的奇异值分解 backsolve：解上三角或下三角方程组 chol：Choleski分解
qr：矩阵的QR分解 chol2inv：由Choleski分解求逆
<，>，<=，>=，==，!=：比较运算符 !，&，&&，|，||，xor()：逻辑运算符 logical：生成逻辑向量 all，any：逻辑向量都为真或存在真
ifelse()：二者择一 match，%in%：查找 unique：找出互不相同的元素 which：找到真值下标集合 duplicated：找到重复元素
optimize，uniroot，polyroot：一维优化与求根
if，else，ifelse，switch：分支 for，while，repeat，break，next：循环 apply，lapply，sapply，tapply，sweep：替代循环的函数。
function：函数定义 source：调用文件 call：函数调用 .C，.Fortran：调用C或者Fortran子程序的动态链接库。 Recall：递归调用
browser，debug，trace，traceback：程序调试 options：指定系统参数 missing：判断虚参是否有对应实参 nargs：参数个数 stop：终止函数执行
on.exit：指定退出时执行 eval，expression：表达式计算 system.time：表达式计算计时 invisible：使变量不显示 menu：选择菜单（字符列表菜单）
其它与函数有关的还有：delay，delete.response，deparse，do.call，dput，environment ，，formals，format.info，interactive，
is.finite，is.function，is.language，is.recursive ，match.arg，match.call，match.fun，model.extract，name，parse，substitute，sys.parent ，warning，machine
cat，print：显示对象 sink：输出转向到指定文件 dump，save，dput，write：输出对象 scan，read.table，load，dget：读入
ls，objects：显示对象列表 rm, remove：删除对象 q，quit：退出系统 .First，.Last：初始运行函数与退出运行函数。
options：系统选项 ?，help，help.start，apropos：帮助功能 data：列出数据集分析
每一种分布有四个函数：d――density（密度函数），p――分布函数，q――分位数函数，r――随机数函数。
比如，正态分布的这四个函数为dnorm，pnorm，qnorm，rnorm。下面我们列出各分布后缀，前面加前缀d、p、q或r就构成函数名：
norm：正态，t：t分布，f：F分布，chisq：卡方（包括非中心） unif：均匀，exp：指数，weibull：威布尔，gamma：伽玛，beta：贝塔
lnorm：对数正态，logis：逻辑分布，cauchy：柯西， binom：二项分布，geom：几何分布，hyper：超几何，nbinom：负二项，pois：泊松 signrank：符号秩，
wilcox：秩和，tukey：学生化极差
sum, mean, var, sd, min, max, range, median, IQR（四分位间距）等为统计量，sort，order，rank与排序有关，其它还有ave，fivenum，mad，quantile，stem等。
R中已实现的有chisq.test，prop.test，t.test。
cor，cov.wt，var：协方差阵及相关阵计算 biplot，biplot.princomp：多元数据biplot图 cancor：典则相关 princomp：主成分分析 hclust：谱系聚类
kmeans：k-均值聚类 cmdscale：经典多维标度 其它有dist，mahalanobis，cov.rob。
ts：时间序列对象 diff：计算差分 time：时间序列的采样时间 window：时间窗
lm，glm，aov：线性模型、广义线性模型、方差
1、RMSE（均方根误差）即标准误差：假如数据在A1：Z1
标准方差用函数=STDEV（A1:Z1）
方差用函数=VARA（A1:Z1）
2、MRE（平均相对误差）
Excel/函数/统计/STDEV（Sd）
计算出标准偏差Sd值，然后除以平均数再×100％就可以了。
为了找到均方根误差，我们首先需要找到残差（也称为误差，我们需要对这些值均方根），然后需要计算这些残差的均方根。因此，如果我们有一个线性回归模型对象说M，则均方根误差可以找到为sqrt（mean（M $residuals ^ 2））。
示例
x1<-rnorm(500,50,5)
y1<-rnorm(500,50,2)
M1<-lm(y1~x1)
summary(M1)
输出结果
Call:
lm(formula = y1 ~ x1)
Residuals:
Min 1QMedian3QMax
-5.6621 -1.2257 -0.0272 1.4151 6.6421
Coefficients:
EstimateStd.Errort value Pr(>|t|)
(Intercept) 50.178943 0.915473 54.812 <2e-16 ***
x1 -0.002153 0.018241 -0.118 0.906
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.966 on 498 degrees of freedom
Multiple R-squared: 2.798e-05, Adjusted R-squared: -0.00198
F-statistic: 0.01393 on 1 and 498 DF, p-value: 0.9061
从模型M1中找到均方根误差-
示例
sqrt(mean(M1$residuals^2))
输出结果
[1] 1.961622
示例
x2<-rnorm(5000,125,21)
y2<-rnorm(5000,137,10)
M2<-lm(y2~x2)
summary(M2)
输出结果
Call:
lm(formula = y2 ~ x2)
Residuals:
Min 1QMedian3QMax
-37.425 -7.005 -0.231 6.836 36.627
Coefficients:
Estimate Std.Error t value Pr(>|t|)
(Intercept) 138.683501 0.851247 162.918 <2e-16 ***
x2 -0.014386 0.006735 -2.136 0.0327 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.06 on 4998 degrees of freedom
Multiple R-squared: 0.0009121, Adjusted R-squared: 0.0007122
F-statistic: 4.563 on 1 and 4998 DF, p-value: 0.03272
从模型M2中找到均方根误差：
示例
sqrt(mean(M2$residuals^2))
输出结果
[1] 10.05584
示例
x37<-rpois(500,5)
y3<-rpois(500,10)
M3<-lm(y3~x3)
summary(M3)
输出结果
Call:
lm(formula = y3 ~ x3)
Residuals:
Min 1QMedian3QMax
-7.9004 -1.9928 -0.2155 2.1921 9.3770
Coefficients:
EstimateStd.Error t value Pr(>|t|)
(Intercept) 10.17770 0.3233031.481<2e-16 ***
x3 -0.09244 0.06145-1.5040.133
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.027 on 498 degrees of freedom
Multiple R-squared: 0.004524, Adjusted R-squared: 0.002525
F-statistic: 2.263 on 1 and 498 DF, p-value: 0.1331
从模型M3查找均方根误差-
示例
sqrt(mean(M3$residuals^2))
输出结果
[1] 3.020734
示例
x4<-runif(50000,5,10)
y4<-runif(50000,2,10)
M4<-lm(y4~x4)
summary(M4)
输出结果
Call:
lm(formula = y4 ~ x4)
Residuals:
Min1Q Median 3QMax
-4.0007 -1.9934 -0.0063 1.9956 3.9995
Coefficients:
EstimateStd.Error t value Pr(>|t|)
(Intercept) 5.9994268 0.0546751 109.729 <2e-16 ***
x40.0001572 0.0071579 0.0220.982
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.309 on 49998 degrees of freedom
Multiple R-squared: 9.646e-09, Adjusted R-squared: -1.999e-05
F-statistic: 0.0004823 on 1 and 49998 DF, p-value: 0.9825
从模型M4找到均方根误差-
示例
sqrt(mean(M4$residuals^2))
输出结果
[1] 2.308586
示例
x5<-sample(5001:9999,100000,replace=TRUE)
y5<-sample(1000:9999,100000,replace=TRUE)
M5<-lm(y5~x5)
summary(M5)
输出结果
Call:
lm(formula = y5 ~ x5)
Residuals:
Min 1QMedian 3Q Max
-4495 -2242-42230 4512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.504e+03 4.342e+01 126.765 <2e-16 ***
x5-1.891e-03 5.688e-03 -0.333 0.74
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2594 on 99998 degrees of freedom
Multiple R-squared: 1.106e-06, Adjusted R-squared: -8.895e-06
F-statistic: 0.1106 on 1 and 99998 DF, p-value: 0.7395
从模型M5中找到均方根误差<
示例
sqrt(mean(M5$residuals^2))
输出结果
[1] 2593.709
