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

计算β多样性指数需要用到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

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

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

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

NMDS

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,723评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,080评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,604评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,440评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,431评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,499评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,893评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,541评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,751评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,547评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,619评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,320评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,890评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,896评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,137评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,796评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,335评论 2 342

推荐阅读更多精彩内容