計(jì)算β多樣性指數(shù)需要用到phyloseq包踏志。它的安裝方式不同于簡(jiǎn)單的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包后羹与,開(kāi)始讀取數(shù)據(jù),前面計(jì)算α多樣性冠摄,用到的是read.table……
qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")
#讀取數(shù)據(jù)梗摇,參數(shù)都是文件名拓哟,注意加后綴
#otufilename指定out表格,mapfilename指定map文件(分組數(shù)據(jù))
#treefilename指定有根進(jìn)化樹(shù)文件
#refseqfilename指定代表序列文件
otu<-qiimedata@otu_table@.Data
#從qiimedata里面提取otu
sum_of_otus<-colSums(t(otu))
#t_轉(zhuǎn)置,colsums計(jì)算列的和,即計(jì)算各個(gè)otu檢測(cè)到的總序列數(shù)伶授,為了篩掉一些總序列數(shù)過(guò)低的otu(可能是測(cè)序錯(cuò)誤)
sum_of_otus
#查看otu總序列數(shù)
selected_otu<-names(sum_of_otus)[sum_of_otus>10]
#獲取總序列數(shù)大于10的otu id
sub_qiimedata <- prune_taxa(selected_otu, qiimedata)
#篩選總序列數(shù)大于10的otu的phyloseq數(shù)據(jù)
weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')
#計(jì)算樣本間加權(quán)unifrac
unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')
#計(jì)算樣本間非加權(quán)unifrac
bray_curtis <- distance(sub_qiimedata, method='bray')
write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)
#保存距離矩陣
#計(jì)算樣本間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排序分析結(jié)果可視化
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加標(biāo)題,stat_ellipse加橢圓
#用scale_colour_manual(values=c())自定義顏色糜烹,可查顏色的16進(jìn)制對(duì)照表
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排序分析結(jié)果可視化
# color=“Group1”指定不同分組的點(diǎn)染不同顏色
p1
p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))
#對(duì)圖片進(jìn)行適當(dāng)修飾违诗, stat_ellipse()加橢圓, ggtitle()加標(biāo)題
ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width
PCoA中的兩個(gè)點(diǎn)距離疮蹦,接近β多樣性指數(shù)
PCA(Principal Components Analysis)即主成分分析诸迟,也稱(chēng)主分量分析或主成分回歸分析法,首先利用線(xiàn)性變換愕乎,將數(shù)據(jù)變換到一個(gè)新的坐標(biāo)系統(tǒng)中;然后再利用降維的思想阵苇,使得任何數(shù)據(jù)投影的第一大方差在第一個(gè)坐標(biāo)(稱(chēng)為第一主成分)上,第二大方差在第二個(gè)坐標(biāo)(第二主成分)上感论。這種降維的思想首先減少數(shù)據(jù)集的維數(shù)绅项,同時(shí)還保持?jǐn)?shù)據(jù)集的對(duì)方差貢獻(xiàn)最大的特征,最終使數(shù)據(jù)直觀呈現(xiàn)在二維坐標(biāo)系比肄。
PCoA(Principal Co-ordinates Analysis)分析即主坐標(biāo)分析快耿,可呈現(xiàn)研究數(shù)據(jù)相似性或差異性的可視化坐標(biāo),是一種非約束性的數(shù)據(jù)降維分析方法薪前,可用來(lái)研究樣本群落組成的相似性或相異性润努。它與PCA類(lèi)似,通過(guò)一系列的特征值和特征向量進(jìn)行排序后示括,選擇主要排在前幾位的特征值铺浇,找到距離矩陣中最主要的坐標(biāo),結(jié)果是數(shù)據(jù)矩陣的一個(gè)旋轉(zhuǎn)垛膝,它沒(méi)有改變樣本點(diǎn)之間的相互位置關(guān)系鳍侣,只是改變了坐標(biāo)系統(tǒng)丁稀。兩者的區(qū)別為PCA是基于樣本的相似系數(shù)矩陣(如歐式距離)來(lái)尋找主成分,而PCoA是基于距離矩陣(歐式距離以外的其他距離)來(lái)尋找主坐標(biāo)倚聚。
NMDS圖中兩個(gè)點(diǎn)的距離的排序线衫,接近β多樣性指數(shù)的排序