R語(yǔ)言計(jì)算β多樣性指數(shù)及分析

計(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

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

NMDS圖中兩個(gè)點(diǎn)的距離的排序线衫,接近β多樣性指數(shù)的排序

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個(gè)濱河市惑折,隨后出現(xiàn)的幾起案子授账,更是在濱河造成了極大的恐慌,老刑警劉巖惨驶,帶你破解...
    沈念sama閱讀 211,639評(píng)論 6 492
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件白热,死亡現(xiàn)場(chǎng)離奇詭異,居然都是意外死亡粗卜,警方通過(guò)查閱死者的電腦和手機(jī)屋确,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,277評(píng)論 3 385
  • 文/潘曉璐 我一進(jìn)店門(mén),熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)续扔,“玉大人攻臀,你說(shuō)我怎么就攤上這事∩疵粒” “怎么了刨啸?”我有些...
    開(kāi)封第一講書(shū)人閱讀 157,221評(píng)論 0 348
  • 文/不壞的土叔 我叫張陵,是天一觀的道長(zhǎng)砌些。 經(jīng)常有香客問(wèn)我呜投,道長(zhǎng),這世上最難降的妖魔是什么存璃? 我笑而不...
    開(kāi)封第一講書(shū)人閱讀 56,474評(píng)論 1 283
  • 正文 為了忘掉前任仑荐,我火速辦了婚禮,結(jié)果婚禮上纵东,老公的妹妹穿的比我還像新娘粘招。我一直安慰自己,他們只是感情好偎球,可當(dāng)我...
    茶點(diǎn)故事閱讀 65,570評(píng)論 6 386
  • 文/花漫 我一把揭開(kāi)白布洒扎。 她就那樣靜靜地躺著,像睡著了一般衰絮。 火紅的嫁衣襯著肌膚如雪袍冷。 梳的紋絲不亂的頭發(fā)上,一...
    開(kāi)封第一講書(shū)人閱讀 49,816評(píng)論 1 290
  • 那天猫牡,我揣著相機(jī)與錄音胡诗,去河邊找鬼。 笑死,一個(gè)胖子當(dāng)著我的面吹牛煌恢,可吹牛的內(nèi)容都是我干的骇陈。 我是一名探鬼主播,決...
    沈念sama閱讀 38,957評(píng)論 3 408
  • 文/蒼蘭香墨 我猛地睜開(kāi)眼瑰抵,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼你雌!你這毒婦竟也來(lái)了?” 一聲冷哼從身側(cè)響起二汛,我...
    開(kāi)封第一講書(shū)人閱讀 37,718評(píng)論 0 266
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤婿崭,失蹤者是張志新(化名)和其女友劉穎,沒(méi)想到半個(gè)月后习贫,有當(dāng)?shù)厝嗽跇?shù)林里發(fā)現(xiàn)了一具尸體逛球,經(jīng)...
    沈念sama閱讀 44,176評(píng)論 1 303
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,511評(píng)論 2 327
  • 正文 我和宋清朗相戀三年苫昌,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片幸海。...
    茶點(diǎn)故事閱讀 38,646評(píng)論 1 340
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡祟身,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出物独,到底是詐尸還是另有隱情袜硫,我是刑警寧澤,帶...
    沈念sama閱讀 34,322評(píng)論 4 330
  • 正文 年R本政府宣布挡篓,位于F島的核電站婉陷,受9級(jí)特大地震影響,放射性物質(zhì)發(fā)生泄漏官研。R本人自食惡果不足惜秽澳,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,934評(píng)論 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望戏羽。 院中可真熱鬧担神,春花似錦、人聲如沸始花。這莊子的主人今日做“春日...
    開(kāi)封第一講書(shū)人閱讀 30,755評(píng)論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)酷宵。三九已至亥贸,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間浇垦,已是汗流浹背炕置。 一陣腳步聲響...
    開(kāi)封第一講書(shū)人閱讀 31,987評(píng)論 1 266
  • 我被黑心中介騙來(lái)泰國(guó)打工, 沒(méi)想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留,地道東北人讹俊。 一個(gè)月前我還...
    沈念sama閱讀 46,358評(píng)論 2 360
  • 正文 我出身青樓垦沉,卻偏偏與公主長(zhǎng)得像,于是被迫代替她去往敵國(guó)和親仍劈。 傳聞我的和親對(duì)象是個(gè)殘疾皇子厕倍,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 43,514評(píng)論 2 348

推薦閱讀更多精彩內(nèi)容