跟著Nature Communications學(xué)數(shù)據(jù)分析:R語言DESeq2包做otu差異豐度分析及火山圖展示結(jié)果

論文

Microbiome differential abundance methods produce different results across 38 datasets

數(shù)據(jù)鏈接

https://figshare.com/articles/dataset/16S_rRNA_Microbiome_Datasets/14531724

代碼鏈接

https://github.com/nearinj/Comparison_of_DA_microbiome_methods

這個人的github主頁還有其他論文的數(shù)據(jù)和代碼

https://github.com/jnmacdonald/differential-abundance-analysis 這個鏈接有很多冠以差異豐度分析的代碼

這兩天在看宏基因組的利用otu豐度數(shù)據(jù)做差異豐度分析,找到了這篇論文,看了摘要,好像是比較了不同差異豐度分析方法獲得結(jié)果的異同。重復(fù)一下這里利用DESeq2做差異豐度分析的代碼

這里我用到的數(shù)據(jù)集是

  • 豐度數(shù)據(jù) ArcticFireSoils_genus_table.tsv
  • 分組數(shù)據(jù) ArcticFireSoils_meta.tsv

這里有一個疑問:論文提供的豐度表格數(shù)據(jù)有兩個歉秫,還有一個是帶rare后綴的,暫時不知道這兩個有啥區(qū)別

首先是讀取數(shù)據(jù)集

ASV_table <- read.table("metagenomics/dat01/ArcticFireSoils_genus_table.tsv", 
                        sep="\t", 
                        skip=1, 
                        header=T, 
                        row.names = 1,
                        comment.char = "", 
                        quote="", check.names = F)
groupings <- read.table("metagenomics/dat01/ArcticFireSoils_meta.tsv", 
                        sep="\t", 
                        row.names = 1, 
                        header=T, 
                        comment.char = "", 
                        quote="", 
                        check.names = F)
dim(ASV_table)
dim(groupings)
groupings$Fire<-factor(groupings$Fire)

這里需要對表示分組的賦予因子,要不然后面deseq2的步驟會有警告信息

判斷一下豐度數(shù)據(jù)的樣本名和分組數(shù)據(jù)的樣本名順序是否一致

identical(colnames(ASV_table), rownames(groupings))

返回false不一致

對兩個樣本名取交集

rows_to_keep <- intersect(colnames(ASV_table), rownames(groupings))

根據(jù)取交集的結(jié)果重新選擇樣本

groupings <- groupings[rows_to_keep,,drop=F]
ASV_table <- ASV_table[,rows_to_keep]

疑問:這里中括號里的drop參數(shù)是啥作用了罢浇?

再次判斷兩個數(shù)據(jù)集中樣本名的順序

identical(colnames(ASV_table), rownames(groupings))

這次返回TRUE

對分組文件的列名進行修改

colnames(groupings)[1] <- "Groupings"

差異豐度分析

library(DESeq2)
dds <- DESeq2::DESeqDataSetFromMatrix(countData = ASV_table,
                                      colData=groupings,
                                      design = ~ Groupings)
dds_res <- DESeq2::DESeq(dds, sfType = "poscounts")

res <- results(dds_res, 
               tidy=T, 
               format="DataFrame",
               contrast = c("Groupings","Fire","Control"))
head(res)
image.png

火山圖代碼

DEG<-res
logFC_cutoff<-2
DEG$change<-as.factor(ifelse(DEG$pvalue<0.05&abs(DEG$log2FoldChange)>logFC_cutoff,
                             ifelse(DEG$log2FoldChange>logFC_cutoff,"UP","DOWN"),
                             "NOT"))
this_title <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                     '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                     '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
DEG<-na.omit(DEG)
library(ggplot2)
ggplot(data=DEG,aes(x=log2FoldChange,
                    y=-log10(pvalue),
                    color=change))+
  geom_point(alpha=0.8,size=3)+
  labs(x="log2 fold change")+ ylab("-log10 pvalue")+
  ggtitle(this_title)+theme_bw(base_size = 20)+
  theme(plot.title = element_text(size=15,hjust=0.5),)+
  scale_color_manual(values=c('#a121f0','#bebebe','#ffad21')) -> p1
p1+xlim(NA,10)+ylim(NA,30) -> p2

library(patchwork)
p1+p2
image.png

今天推文的示例數(shù)據(jù)和代碼可以在公眾號查看獲取方式

歡迎大家關(guān)注我的公眾號

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號 主要分享:1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子沐祷;2嚷闭、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)赖临、群體遺傳學(xué)文獻閱讀筆記胞锰;3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記兢榨!

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末嗅榕,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子吵聪,更是在濱河造成了極大的恐慌凌那,老刑警劉巖,帶你破解...
    沈念sama閱讀 223,002評論 6 519
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件吟逝,死亡現(xiàn)場離奇詭異帽蝶,居然都是意外死亡,警方通過查閱死者的電腦和手機块攒,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,357評論 3 400
  • 文/潘曉璐 我一進店門励稳,熙熙樓的掌柜王于貴愁眉苦臉地迎上來佃乘,“玉大人,你說我怎么就攤上這事驹尼∪け埽” “怎么了?”我有些...
    開封第一講書人閱讀 169,787評論 0 365
  • 文/不壞的土叔 我叫張陵扶欣,是天一觀的道長鹅巍。 經(jīng)常有香客問我,道長料祠,這世上最難降的妖魔是什么骆捧? 我笑而不...
    開封第一講書人閱讀 60,237評論 1 300
  • 正文 為了忘掉前任,我火速辦了婚禮髓绽,結(jié)果婚禮上敛苇,老公的妹妹穿的比我還像新娘。我一直安慰自己顺呕,他們只是感情好枫攀,可當(dāng)我...
    茶點故事閱讀 69,237評論 6 398
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著株茶,像睡著了一般来涨。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上启盛,一...
    開封第一講書人閱讀 52,821評論 1 314
  • 那天蹦掐,我揣著相機與錄音,去河邊找鬼僵闯。 笑死卧抗,一個胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的鳖粟。 我是一名探鬼主播社裆,決...
    沈念sama閱讀 41,236評論 3 424
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼向图!你這毒婦竟也來了泳秀?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 40,196評論 0 277
  • 序言:老撾萬榮一對情侶失蹤榄攀,失蹤者是張志新(化名)和其女友劉穎晶默,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體航攒,經(jīng)...
    沈念sama閱讀 46,716評論 1 320
  • 正文 獨居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 38,794評論 3 343
  • 正文 我和宋清朗相戀三年趴梢,在試婚紗的時候發(fā)現(xiàn)自己被綠了漠畜。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片币他。...
    茶點故事閱讀 40,928評論 1 353
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖憔狞,靈堂內(nèi)的尸體忽然破棺而出蝴悉,到底是詐尸還是另有隱情,我是刑警寧澤瘾敢,帶...
    沈念sama閱讀 36,583評論 5 351
  • 正文 年R本政府宣布拍冠,位于F島的核電站,受9級特大地震影響簇抵,放射性物質(zhì)發(fā)生泄漏庆杜。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 42,264評論 3 336
  • 文/蒙蒙 一碟摆、第九天 我趴在偏房一處隱蔽的房頂上張望晃财。 院中可真熱鬧,春花似錦典蜕、人聲如沸断盛。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,755評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽钢猛。三九已至,卻和暖如春轩缤,著一層夾襖步出監(jiān)牢的瞬間命迈,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,869評論 1 274
  • 我被黑心中介騙來泰國打工典奉, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留躺翻,地道東北人。 一個月前我還...
    沈念sama閱讀 49,378評論 3 379
  • 正文 我出身青樓卫玖,卻偏偏與公主長得像公你,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子假瞬,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 45,937評論 2 361

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