TCGA 02 不同數(shù)據(jù)DEG結(jié)果取交集火山圖及韋恩圖展示

1. 三個(gè)數(shù)據(jù)集差異基因火山圖（每10個(gè)腫瘤和正常組織樣本）

導(dǎo)入表達(dá)矩陣，選擇NT組的前10和TP組的前10個(gè)樣本進(jìn)行差異分析

library(TCGAbiolinks)
# 導(dǎo)入dataFilt表達(dá)矩陣
load("dataFilt.RData")

# selection of normal samples "NT"
samplesNT <- TCGAquery_SampleTypes(barcode = colnames(dataFilt),
                                   typesample = c("NT"))
# selection of tumor samples "TP"
samplesTP <- TCGAquery_SampleTypes(barcode = colnames(dataFilt), 
                                   typesample = c("TP"))
# Diff.expr.analysis (DEA)
DEG.LUAD.edgeR <- TCGAanalyze_DEA(mat1 = dataFilt[,samplesNT[1:10]],
                                  mat2 = dataFilt[,samplesTP[1:10]],
                                  pipeline="edgeR",
                                  batch.factors = c("TSS"),
                                  Cond1type = "Normal",
                                  Cond2type = "Tumor",
                                  voom =FALSE, 
                                  method = "glmLRT",
                                  # fdr.cut = 0.01,  #保留FDR<0.01的基因
                                  # logFC.cut = 1 #保留logFC＞1的基因
)
# ----------------------- DEA -------------------------------
#   there are Cond1 type Normal in  10 samples
# there are Cond2 type Tumor in  10 samples
# there are  12980 features as miRNA or genes 
# I Need about  8.7 seconds for this DEA. [Processing 30k elements /s]  
# ----------------------- END DEA -------------------------------

繪制第一個(gè)火山圖

valcano_data <- data.frame(genes=rownames(DEG.LUAD.edgeR), 
                           logFC=DEG.LUAD.edgeR$logFC, 
                           FDR=DEG.LUAD.edgeR$FDR,
                           group=rep("NotSignificant", 
                                     nrow(DEG.LUAD.edgeR)),
                           stringsAsFactors = F)

valcano_data[which(valcano_data['FDR'] < 0.05 & 
                     valcano_data['logFC'] > 1.5),"group"] <- "Increased"
valcano_data[which(valcano_data['FDR'] < 0.05 &
                     valcano_data['logFC'] < -1.5),"group"] <- "Decreased"

cols = c("darkgrey","#00B2FF","orange")
names(cols) = c("NotSignificant","Increased","Decreased")

library(ggplot2)
vol1 <- ggplot(valcano_data, aes(x = logFC, y = -log10(FDR), color = group))+
  scale_colour_manual(values = cols) +
  ggtitle(label = "Volcano Plot 1", subtitle = "LUAD 1-10 samples volcano plot") +
  geom_point(size = 2.5, alpha = 1, na.rm = T) +
  theme_bw(base_size = 14) + 
  theme(legend.position = "right") + 
  xlab(expression(log[2]("logFC"))) + 
  ylab(expression(-log[10]("FDR"))) +
  geom_hline(yintercept = 1.30102, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = 1.5849, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = -1.5849, colour="#990000", linetype="dashed")+ 
  scale_y_continuous(trans = "log1p")

選擇NT組的11-20和TP組的11-20樣本進(jìn)行差異分析

# Diff.expr.analysis (DEA)
DEG.LUAD.edgeR2 <- TCGAanalyze_DEA(mat1 = dataFilt[,samplesNT[11:20]],
                                  mat2 = dataFilt[,samplesTP[11:20]],
                                  pipeline="edgeR",
                                  batch.factors = c("TSS"),
                                  Cond1type = "Normal",
                                  Cond2type = "Tumor",
                                  voom =FALSE, 
                                  method = "glmLRT",
                                  # fdr.cut = 0.01,  #保留FDR<0.01的基因
                                  # logFC.cut = 1 #保留logFC＞1的基因
)
# ----------------------- DEA -------------------------------
#   there are Cond1 type Normal in  10 samples
# there are Cond2 type Tumor in  10 samples
# there are  12980 features as miRNA or genes 
# I Need about  8.7 seconds for this DEA. [Processing 30k elements /s]  
# ----------------------- END DEA -------------------------------

繪制第2個(gè)火山圖

valcano_data2 <- data.frame(genes=rownames(DEG.LUAD.edgeR2), 
                           logFC=DEG.LUAD.edgeR2$logFC, 
                           FDR=DEG.LUAD.edgeR2$FDR,
                           group=rep("NotSignificant", 
                                     nrow(DEG.LUAD.edgeR2)),
                           stringsAsFactors = F)

valcano_data2[which(valcano_data2['FDR'] < 0.05 & 
                     valcano_data2['logFC'] > 1.5),"group"] <- "Increased"
valcano_data2[which(valcano_data2['FDR'] < 0.05 &
                     valcano_data2['logFC'] < -1.5),"group"] <- "Decreased"

cols = c("darkgrey","#00B2FF","orange")
names(cols) = c("NotSignificant","Increased","Decreased")

library(ggplot2)
vol2 <- ggplot(valcano_data2, aes(x = logFC, y = -log10(FDR), color = group))+
  scale_colour_manual(values = cols) +
  ggtitle(label = "Volcano Plot 2", subtitle = "LUAD 11-20 volcano plot") +
  geom_point(size = 2.5, alpha = 1, na.rm = T) +
  theme_bw(base_size = 14) + 
  theme(legend.position = "right") + 
  xlab(expression(log[2]("logFC"))) + 
  ylab(expression(-log[10]("FDR"))) +
  geom_hline(yintercept = 1.30102, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = 1.5849, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = -1.5849, colour="#990000", linetype="dashed")+ 
  scale_y_continuous(trans = "log1p")

選擇NT組的21-30和TP組的21-30樣本進(jìn)行差異分析

# Diff.expr.analysis (DEA)
DEG.LUAD.edgeR3 <- TCGAanalyze_DEA(mat1 = dataFilt[,samplesNT[21:30]],
                                  mat2 = dataFilt[,samplesTP[21:30]],
                                  pipeline="edgeR",
                                  batch.factors = c("TSS"),
                                  Cond1type = "Normal",
                                  Cond2type = "Tumor",
                                  voom =FALSE, 
                                  method = "glmLRT",
                                  # fdr.cut = 0.01,  #保留FDR<0.01的基因
                                  # logFC.cut = 1 #保留logFC＞1的基因
)
# ----------------------- DEA -------------------------------
#   there are Cond1 type Normal in  10 samples
# there are Cond2 type Tumor in  10 samples
# there are  12980 features as miRNA or genes 
# I Need about  8.7 seconds for this DEA. [Processing 30k elements /s]  
# ----------------------- END DEA -------------------------------

繪制第3個(gè)火山圖

valcano_data3 <- data.frame(genes=rownames(DEG.LUAD.edgeR3), 
                            logFC=DEG.LUAD.edgeR3$logFC, 
                            FDR=DEG.LUAD.edgeR3$FDR,
                            group=rep("NotSignificant", 
                                      nrow(DEG.LUAD.edgeR3)),
                            stringsAsFactors = F)

valcano_data3[which(valcano_data3['FDR'] < 0.05 & 
                     valcano_data3['logFC'] > 1.5),"group"] <- "Increased"
valcano_data3[which(valcano_data3['FDR'] < 0.05 &
                     valcano_data3['logFC'] < -1.5),"group"] <- "Decreased"

cols = c("darkgrey","#00B2FF","orange")
names(cols) = c("NotSignificant","Increased","Decreased")

library(ggplot2)
vol3 <- ggplot(valcano_data3, aes(x = logFC, y = -log10(FDR), color = group))+
  scale_colour_manual(values = cols) +
  ggtitle(label = "Volcano Plot 3", subtitle = "LUAD 21-30 volcano plot") +
  geom_point(size = 2.5, alpha = 1, na.rm = T) +
  theme_bw(base_size = 14) + 
  theme(legend.position = "right") + 
  xlab(expression(log[2]("logFC"))) + 
  ylab(expression(-log[10]("FDR"))) +
  geom_hline(yintercept = 1.30102, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = 1.5849, colour="#990000", linetype="dashed") + 
  geom_vline(xintercept = -1.5849, colour="#990000", linetype="dashed")+ 
  scale_y_continuous(trans = "log1p")
library(cowplot)
library(patchwork)
vol1+vol2+vol3
3個(gè)火山圖合并

2. upset圖和韋恩圖分析

提取3個(gè)數(shù)據(jù)的差異基因列表

DEG1_up <- valcano_data[valcano_data$group=="Increased", "genes"]
DEG1_down <- valcano_data[valcano_data$group=="Decreased", "genes"]

DEG2_up <- valcano_data[valcano_data2$group=="Increased", "genes"]
DEG2_down <- valcano_data[valcano_data2$group=="Decreased", "genes"]

DEG3_up <- valcano_data[valcano_data3$group=="Increased", "genes"]
DEG3_down <- valcano_data[valcano_data3$group=="Decreased", "genes"]

用Y叔開發(fā)的ggupset做Upset圖

# devtools::install_github("GuangchuangYu/yyplot")
# devtools::install_github("GuangchuangYu/UpSetR")
# install.packages("venneuler")
# remove.packages("ggplot2")
# install.packages("ggimage")
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("ComplexHeatmap")

# 上調(diào)基因
lt_up = list(TCGA_1 = DEG1_up,
             TCGA_2 = DEG2_up,
             TCGA_3 = DEG3_up)
dat_up<- ComplexHeatmap::list_to_matrix(lt_up)
dat_plot_up <-  data.frame(dat_up)

require(UpSetR)
p1 <- upset(dat_plot_up, sets=c("TCGA_3", "TCGA_2", "TCGA_1"),
  keep.order = TRUE)
require(ggplotify)
g1 <- as.ggplot(p1) + ggtitle("Up regulated")


# 下調(diào)基因
lt_down = list(TCGA_1 = DEG1_down,
               TCGA_2 = DEG2_down,
               TCGA_3 = DEG3_down)
dat_down<- ComplexHeatmap::list_to_matrix(lt_down)
dat_plot_down <-  data.frame(dat_down)

require(UpSetR)
p2 <- upset(dat_plot_down, , sets=c("TCGA_3", "TCGA_2", "TCGA_1"),
            keep.order = TRUE)
require(ggplotify)
g2 <- as.ggplot(p2) + ggtitle("Down regulated")

# 拼圖
library(cowplot)
library(patchwork)
g1+g2
upset圖合并

做韋恩圖

require(yyplot)
library("ggsci")
g3 <- ggvenn(dat_plot_up) + theme_void() +  
  scale_fill_jco() + ggtitle("Up regulated genes")
g4 <- ggvenn(dat_plot_down) + theme_void() +  
  scale_fill_jco() + ggtitle("Down regulated genes")
g3 + g4
韋恩圖合并

upset和韋恩圖拼接

require(ggimage)
g1_3 <- g1 + geom_subview(subview=g3+theme_void(), x=.78, y=.8, w=.5, h=.5)
g2_4 <- g2 + geom_subview(subview=g4+theme_void(), x=.78, y=.8, w=.5, h=.5)

g1_3+g2_4
upset和韋恩圖合并

3. 火山圖、upset圖+韋恩圖合并

(vol1|vol2|vol3)/
  (g1_3|g2_4)
火山圖-韋恩圖-upset圖合并

挑選出三個(gè)數(shù)據(jù)集中共同上調(diào)或下調(diào)的基因

up_common <- Reduce(intersect, list(DEG1_up, DEG2_up, DEG3_up))
down_common <- Reduce(intersect, list(DEG1_down, DEG2_down, DEG3_down))
background_genes <- Reduce(union, list(DEG1_up, DEG2_up, DEG3_up,
                                       DEG1_down, DEG2_down, DEG3_down))
up_common_df <- data.frame(gene=up_common, 
                           logFC=valcano_data[valcano_data$genes %in% 
                                                up_common, 
                                              "logFC"])
down_common_df <- data.frame(gene=down_common, 
                             logFC=valcano_data[valcano_data$genes %in% 
                                                  down_common, 
                                                "logFC"])
background_genes_df <- data.frame(gene=background_genes, 
                              logFC=valcano_data[valcano_data$genes %in% 
                                                   background_genes, 
                                                 "logFC"], stringsAsFactors = F)
all_gene_df <- valcano_data[, c("genes", "logFC")]
save(list=c("up_common_df", "down_common_df", "background_genes_df", "all_gene_df"), 
     file="filtered_genes.RData")
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌，老刑警劉巖鹊碍，帶你破解...
    沈念sama閱讀 206,013評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異食绿，居然都是意外死亡侈咕，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門器紧，熙熙樓的掌柜王于貴愁眉苦臉地迎上來耀销，“玉大人，你說我怎么就攤上這事铲汪⌒芪荆” “怎么了罐柳？”我有些...
    開封第一講書人閱讀 152,370評論 0贊 342
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵膏萧，是天一觀的道長息尺。 經(jīng)常有香客問我，道長脚囊，這世上最難降的妖魔是什么转晰？ 我笑而不...
    開封第一講書人閱讀 55,168評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮士飒，結(jié)果婚禮上查邢，老公的妹妹穿的比我還像新娘。我一直安慰自己酵幕，他們只是感情好扰藕，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,153評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著芳撒，像睡著了一般邓深。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上笔刹，一...
    開封第一講書人閱讀 48,954評論 1贊 283
  • 城市分裂傳說
    那天芥备，我揣著相機(jī)與錄音，去河邊找鬼舌菜。 笑死萌壳，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的日月。 我是一名探鬼主播袱瓮，決...
    沈念sama閱讀 38,271評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼爱咬！你這毒婦竟也來了尺借？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 36,916評論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤精拟，失蹤者是張志新（化名）和其女友劉穎燎斩，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體串前，經(jīng)...
    沈念sama閱讀 43,382評論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡瘫里，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,877評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了荡碾。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片谨读。...
    茶點(diǎn)故事閱讀 37,989評論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡，死狀恐怖坛吁，靈堂內(nèi)的尸體忽然破棺而出劳殖，到底是詐尸還是另有隱情铐尚，我是刑警寧澤，帶...
    沈念sama閱讀 33,624評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布哆姻，位于F島的核電站宣增，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏矛缨。R本人自食惡果不足惜爹脾，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,209評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望箕昭。 院中可真熱鬧灵妨，春花似錦、人聲如沸落竹。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0贊 19
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽述召。三九已至朱转，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間积暖，已是汗流浹背藤为。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留呀酸，地道東北人凉蜂。 一個(gè)月前我還...
    沈念sama閱讀 45,401評論 2贊 352
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像性誉，于是被迫代替她去往敵國和親窿吩。 傳聞我的和親對象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,700評論 2贊 345