install.packages("ggplot2")

install.packages("ggpubr")

install.packages("vegan")

计算Shannon-香农指数和Simpson-辛普森指数的命令在vegan包中，计算各组显著性的命令在ggpubr包中；画图使用ggplot命令，在行使每个命令之前一定要加载相应的包，如下:

library(ggplot2)

library(ggpubr)

library(vegan)

拿到一个otu表格，要先计算香农指数和辛普森指数，操作如下:

otu=read.table('D:/r-working/feature-table.taxonomy.txt',row.names = 1,skip=1,header=T,comment.char ='',sep='\t')

#读取out表格

#'D:/feature table.taxonomy.txt'为文件路径，注意斜线方向

#row.names = 1指定第一列为行名

#skip=1跳过第一行不读

#header=T指定第一个有效行为列名

#sep='\t'表示指定制表符为分隔符

#comment.char=''表示设置注释符号为空字符‘’，这样#后面的内容就不会被省略

otu=otu[,-ncol(otu)]

#去除表格的最后一列，无用信息

otu=t(otu)

#表格转置，必须将样品名作为行名

shannon=diversity(otu,"shannon")

#计算香农指数,先加载vegan包

shannon

#查看香农指数

simpson=diversity(otu,"simpson")

#计算辛普森指数,先加载vegan包

simpson

#查看辛普森指数

alpha=data.frame(shannon,simpson,check.names=T)

#合并两个指数

write.table(alpha,"D:/r-working/alpha-summary.xls",sep='\t',quote=F)

#存储数据，注意路径使用反斜杠

将各样本进行分组，并进行画图，操作如下:

map<-read.table('D:/r-working/mapping_file.txt',row.names = 1,header=T,comment.char ='',sep='\t',check.names=F)

#读取分组表格

group<-map["Group1"]

#提取需要的分组,'Group1'是表中的分组列名，包括A,B,C三组

alpha<-alpha[match(rownames(group),rownames(alpha)),]

#重排alpha的行的顺序,使其与group的样本id(行名)一致

data<-data.frame(group,alpha,check.rows=T)

#合并两个表格.'<-'与'='同属赋值的含义.

p=ggplot(data=data,aes(x=Group1,y=shannon))+geom_boxplot(fill=rainbow(7)[2])

#data = data指定数据表格

#x=Group1指定作为x轴的数据列名

#y=shannon指定作为y轴的数据列名

#geom_boxplot()表示画箱线图

#fill=rainbow(7)[2]指定填充色

此处用到ggplot2包画箱线图，将画图函数赋值给p后，可以用‘+’不断进行图层叠加，给图片p增加新的特性

p

#查看p

mycompare=list(c('A','B'),c('A','C'),c('B','C'))

#指定多重比较的分组对

mycompare

p<-p+stat_compare_means(comparisons=mycompare,label = "p.signif",method = 'wilcox')

#添加显著性标记的第一种方法，在此之前先加载ggpubr包

p<-p+ylim(2,5.5)

#调整图像的外观

本文基于R语言进行基本数据统计分析，包括基本作图，线性拟合，逻辑回归，bootstrap采样和Anova方差分析的实现及应用。

不多说，直接上代码，代码中有注释。

1. 基本作图（盒图，qq图）

#basic plot

boxplot(x)

qqplot(x,y)

2. 线性拟合

#linear regression

n = 10

x1 = rnorm(n)#variable 1

x2 = rnorm(n)#variable 2

y = rnorm(n)*3

mod = lm(y~x1+x2)

model.matrix(mod) #erect the matrix of mod

plot(mod) #plot residual and fitted of the solution, Q-Q plot and cook distance

summary(mod) #get the statistic information of the model

hatvalues(mod) #very important, for abnormal sample detection

3. 逻辑回归

#logistic regression

x <- c(0, 1, 2, 3, 4, 5)

y <- c(0, 9, 21, 47, 60, 63) # the number of successes

n <- 70 #the number of trails

z <- n - y #the number of failures

b <- cbind(y, z) # column bind

fitx <- glm(b~x,family = binomial) # a particular type of generalized linear model

print(fitx)

plot(x,y,xlim=c(0,5),ylim=c(0,65)) #plot the points (x,y)

beta0 <- fitx$coef[1]

beta1 <- fitx$coef[2]

fn <- function(x) n*exp(beta0+beta1*x)/(1+exp(beta0+beta1*x))

par(new=T)

curve(fn,0,5,ylim=c(0,60)) # plot the logistic regression curve

3. Bootstrap采样

# bootstrap

# Application: 随机采样，获取最大eigenvalue占所有eigenvalue和之比，并画图显示distribution

dat = matrix(rnorm(100*5),100,5)

no.samples = 200 #sample 200 times

# theta = matrix(rep(0,no.samples*5),no.samples,5)

theta =rep(0,no.samples*5)

for (i in 1:no.samples)

{

j = sample(1:100,100,replace = TRUE)#get 100 samples each time

datrnd = dat[j,]#select one row each time

lambda = princomp(datrnd)$sdev^2#get eigenvalues

# theta[i,] = lambda

theta[i] = lambda[1]/sum(lambda)#plot the ratio of the biggest eigenvalue

}

# hist(theta[1,]) #plot the histogram of the first(biggest) eigenvalue

hist(theta)#plot the percentage distribution of the biggest eigenvalue

sd(theta)#standard deviation of theta

#上面注释掉的语句，可以全部去掉注释并将其下一条语句注释掉，完成画最大eigenvalue分布的功能

4. ANOVA方差分析

#Application：判断一个自变量是否有影响 (假设我们喂3种维他命给3头猪，想看喂维他命有没有用)

#

y = rnorm(9)#weight gain by pig(Yij, i is the treatment, j is the pig_id), 一般由用户自行输入

#y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)

Treatment <- factor(c(1,2,3,1,2,3,1,2,3)) #each {1,2,3} is a group

mod = lm(y~Treatment) #linear regression

print(anova(mod))

#解释：Df（degree of freedom）

#Sum Sq: deviance (within groups, and residuals) 总偏差和

# Mean Sq: variance (within groups, and residuals) 平均方差和

# compare the contribution given by Treatment and Residual

#F value: Mean Sq(Treatment)/Mean Sq(Residuals)

#Pr(>F): p-value. 根据p-value决定是否接受Hypothesis H0：多个样本总体均数相等(检验水准为0.05)

qqnorm(mod$residual) #plot the residual approximated by mod

#如果qqnorm of residual像一条直线，说明residual符合正态分布，也就是说Treatment带来的contribution很小，也就是说Treatment无法带来收益（多喂维他命少喂维他命没区别）

如下面两图分别是

（左）用 y = matrix(c(1,10,1,2,10,2,1,9,1),9,1)和

（右）y = rnorm(9)

的结果。可见如果给定猪吃维他命2后体重特别突出的数据结果后，qq图种residual不在是一条直线，换句话说residual不再符合正态分布，i.e., 维他命对猪的体重有影响。