通过3个数据集取差异基因的交集复现一篇文中的韦恩图

老大让复现的是一张韦恩图,来自3个数据集,分别找每个数据集的差异基因,然后取UP和DOWN的交集。原文链接:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6054764/#.

原文中的图片如下

image-20191121212758565

想先用循环下载含有表达矩阵和临床信息的一个list,然后试着用循环

library(GEOquery)

ls<-c('GSE38959','GSE45827','GSE65194')

for (i in 1:length(ls)){

gset[[i]] <- getGEO(ls[i], #系列编号

destdir = '.', #当前目录

getGPL = F)

}

1

2

3

4

5

6

7

但是下载下来的就只是表达矩阵,并不是含有临床信息的list,下载后的表达矩阵如下

image-20191121212900565

接着想再尝试循环读取上面的表达矩阵

tmp<-dir(pattern = ".gz")

tmp

dat<-list()

for (i in 1:length(tmp)) {

dat[[i]]<-read.table(tmp[i],comment.char = '!',sep = '\t',header = T,row.names = 1)

}

1

2

3

4

5

6

上面就能得到含有三个data.frame的list了,上面两个循环就只会到这里,因为还要再去下载数据集的临床信息,所以就分三个数据集分别下载,得到上下调基因然后做韦恩图。

image-20191121212816589

  • 获得第一个数据集的差异基因,下面是代码

下载数据

rm(list = ls())

options(stringsAsFactors = F)

f='GSE38959_eSet.Rdata'

# https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38959

library(GEOquery)

if(!file.exists(f)){

gset <- getGEO('GSE38959', destdir=".",

AnnotGPL = F,

getGPL = F)

save(gset,file=f)

}

load('GSE38959_eSet.Rdata')

class(gset)

length(gset)

class(gset[[1]])

a=gset[[1]]

dat=exprs(a) #a现在是一个对象,取a这个对象通过看说明书知道要用exprs这个函数

dim(dat)

dat[1:4,1:4]

pd=pData(a) #通过查看说明书知道取对象a里的临床信息用pData

tnbc=rownames(pd[grepl('triple',as.character(pd$title)),])

normal=rownames(pd[grepl('normal',as.character(pd$title)),])

dat=dat[,c(tnbc,normal)]

group_list=c(rep('tnbc',length(tnbc)),

rep('normal',length(normal)))

table(group_list)

dat[1:4,1:4]

if(F){

# library(GEOquery)

# #Download GPL file, put it in the current directory, and load it:

# gpl <- getGEO('GPL4133', destdir=".")

# colnames(Table(gpl))

# head(Table(gpl)[,c(1,10)]) ## you need to check this , which column do you need

# probe2gene=Table(gpl)[,c(1,10)]

b=read.table('GPL4133-12599.txt',

sep = '\t',quote = '',

fill=T,header = T)

colnames(b)

probe2gene=b[,c(1,10)]

probe2gene=probe2gene[probe2gene$GENE_SYMBOL != '',]

head(probe2gene)

save(probe2gene,file='probe2gene.Rdata')

}

load(file='probe2gene.Rdata')

head(probe2gene)

ids=probe2gene

colnames(ids)=c('probe_id','symbol')

ids=ids[ids$probe_id %in% rownames(dat),]

ids$probe_id=as.character(ids$probe_id)

dat[1:4,1:4]

dat=dat[ids$probe_id,]

ids$median=apply(dat,1,median) #ids新建median这一列,列名为median,同时对dat这个矩阵按行操作,取每一行的中位数,将结果给到median这一列的每一行

ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序,将对应的行赋值为一个新的ids

ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项,'!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果s

dat=dat[ids$probe_id,] #新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat

rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名

dat[1:4,1:4] #保留每个基因ID第一次出现的信息

dim(dat)

boxplot(dat)

dat=log(dat+1)

dat[1:4,1:4]

save(dat,group_list,file = 'step1-output.Rdata')

load(file = 'step1-output.Rdata')

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

检查数据

rm(list = ls())

options(stringsAsFactors = F)

load(file = 'step1-output.Rdata')

dat[1:4,1:4]

## 下面是画PCA的必须操作,需要看说明书。

dat=t(dat)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换

dat=as.data.frame(dat)

dat=cbind(dat,group_list)

library("FactoMineR")

library("factoextra")

dat.pca <- PCA(dat[,ncol(dat)], graph = FALSE)#现在dat最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的

fviz_pca_ind(dat.pca,

geom.ind = "point",

col.ind = dat$group_list,

palette = c("#00AFBB", "#E7B800"),

addEllipses = TRUE,

legend.title = "Groups"

)

ggsave('all_samples_PCA_by_pCR.png')

rm(list = ls())

load(file = 'step1-output.Rdata')

dat[1:4,1:4]

cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行('1'是按行取,'2'是按列取)取每一行的方差,从小到大排序,取最大的1000个

library(pheatmap)

pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对那些提取出来的1000个基因所在的每一行取出,组合起来为一个新的表达矩阵

n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化

n[n>2]=2

n[n< 2]= 2

n[1:4,1:4]

pheatmap(n,show_colnames =F,show_rownames = F)

ac=data.frame(g=group_list)

rownames(ac)=colnames(n)

pheatmap(n,show_colnames =F,show_rownames = F,

annotation_col=ac,filename = 'heatmap_top1000_sd.png')

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

差异分析

rm(list = ls())

options(stringsAsFactors = F)

load(file = 'step1-output.Rdata')

dat[1:4,1:4]

table(group_list)

boxplot(dat[1,]~group_list)

bp=function(g){

library(ggpubr)

df=data.frame(gene=g,stage=group_list)

p <- ggboxplot(df, x = "stage", y = "gene",

color = "stage", palette = "jco",

add = "jitter")

# Add p-value

p + stat_compare_means()

}

bp(dat[1,]) ## 调用上面定义好的函数,避免同样的绘图代码重复多次敲。

bp(dat[2,])

bp(dat[3,])

bp(dat[4,])

dim(dat)

library(limma)

design=model.matrix(~factor( group_list ))

fit=lmFit(dat,design)

fit=eBayes(fit)

options(digits = 4)

topTable(fit,coef=2,adjust='BH')

## 但是上面的用法做不到随心所欲的指定任意两组进行比较

design <- model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

head(design)

exprSet=dat

rownames(design)=colnames(exprSet)

design

contrast.matrix<-makeContrasts("tnbc-normal",

levels = design)

contrast.matrix ##这个矩阵声明,我们要把 Tumor 组跟 Normal 进行差异分析比较

colnames(design)

deg = function(exprSet,design,contrast.matrix){

fit <- lmFit(exprSet,design)

fit2 <- contrasts.fit(fit, contrast.matrix)

fit2 <- eBayes(fit2)

tempOutput = topTable(fit2, coef=1, n=Inf)

nrDEG = na.omit(tempOutput)

#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)

head(nrDEG)

return(nrDEG)

}

deg = deg(exprSet,design,contrast.matrix)

head(deg)

save(deg,file = 'deg.Rdata')

load(file = 'deg.Rdata')

head(deg)

bp(dat[rownames(deg)[1],])

## for volcano

if(T){

nrDEG=deg

head(nrDEG)

attach(nrDEG)

plot(logFC,log10(P.Value))

library(ggpubr)

df=nrDEG

df$v= log10(P.Value) #df新增加一列'v',值为-log10(P.Value)

ggscatter(df, x = "logFC", y = "v",size=0.5)

df$g=ifelse(df$P.Value>0.01,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因

ifelse( df$logFC >1,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因

ifelse( df$logFC < 1,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1.5 的基因,再if 判断:如果logFC <1.5,则为down(下调)基因,否则为stable基因

)

table(df$g)

df$name=rownames(df)

head(df)

ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')

ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,

label = "name", repel = T,

#label.select = rownames(df)[df$g != 'stable'] ,

label.select =head(rownames(deg)), #挑选一些基因在图中显示出来

palette = c("#00AFBB", "#E7B800", "#FC4E07") )

ggsave('volcano.png')

ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)

df$p_c = ifelse(df$P.Value<0.001,'p<0.001',

ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))

table(df$p_c )

ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,

palette = c("green", "red", "black") )

ggsave('MA.png')

}

## for heatmap

if(T){

load(file = 'step1-output.Rdata')

# 每次都要检测数据

dat[1:4,1:4]

table(group_list)

x=deg$logFC #deg取logFC这列并将其重新赋值给x

names(x)=rownames(deg) #deg取probe_id这列,并将其作为名字给x

cg=c(names(head(sort(x),100)),#对x进行从小到大排列,取前100及后100,并取其对应的探针名,作为向量赋值给cg

names(tail(sort(x),100)))

library(pheatmap)

pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行,所得到的矩阵来画热图

n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于scale这个函数应用在不同组数据间存在差异时,需要行名为样本,因此需要用t(dat[cg,])来转换,最后再转换回来

n[n>2]=2

n[n< 2]= 2

n[1:4,1:4]

pheatmap(n,show_colnames =F,show_rownames = F)

ac=data.frame(g=group_list)

rownames(ac)=colnames(n) #将ac的行名也就分组信息(是‘no TNBC’还是‘TNBC’)给到n的列名,即热图中位于上方的分组信息

pheatmap(n,show_colnames =F,

show_rownames = F,

cluster_cols = T,

annotation_col=ac,filename = 'heatmap_top200_DEG.png') #列名注释信息为ac即分组信息

}

write.csv(deg,file = 'deg.csv')

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

上面都是获得了每个表达矩阵的差异分析结果,上面的代码在另外两个数据集再走一波,当然表达矩阵和分组信息需要自己制作好。每一个数据集得到的差异分析的“deg.Rdata”都保存在自己的文件夹,即使都叫“deg.Rdata”也没关系。下面就是做韦恩图,获得3个数据集的UP及DOWN基因的交集,筛选阈值可以按照文章中或自己调整

  • 韦恩图的代码如下

rm(list = ls()) ## 魔幻操作,一键清空~

# P<0.05 and |logFC|≥2.0,

# GSE38959, 515 upregulated genes and 337 downregulated genes.

# GSE45827, 2,117 genes were upregulated, and 878 genes were downregulated.

# GSE65194, 2,130 upregulated genes and 901 downregulated genes were identified.

load(file = '../GSE38959/deg.Rdata')

head(deg)

## 不同的阈值,筛选到的差异基因数量就不一样,后面的超几何分布检验结果就大相径庭。

logFC_t=2

deg$g=ifelse(deg$P.Value>0.05,'stable',

ifelse( deg$logFC > logFC_t,'UP',

ifelse( deg$logFC < logFC_t,'DOWN','stable') )

)

table(deg$g)

up_GSE38959=rownames(deg[deg$g=='UP',])

down_GSE38959=rownames(deg[deg$g=='DOWN',])

load(file = '../GSE45827/deg.Rdata')

head(deg)

## 不同的阈值,筛选到的差异基因数量就不一样,后面的超几何分布检验结果就大相径庭。

logFC_t=2

deg$g=ifelse(deg$P.Value>0.05,'stable',

ifelse( deg$logFC > logFC_t,'UP',

ifelse( deg$logFC < logFC_t,'DOWN','stable') )

)

table(deg$g)

up_GSE45827=rownames(deg[deg$g=='UP',])

down_GSE45827=rownames(deg[deg$g=='DOWN',])

load(file = '../GSE65194/deg.Rdata')

head(deg)

## 不同的阈值,筛选到的差异基因数量就不一样,后面的超几何分布检验结果就大相径庭。

logFC_t=2

deg$g=ifelse(deg$P.Value>0.05,'stable',

ifelse( deg$logFC > logFC_t,'UP',

ifelse( deg$logFC < logFC_t,'DOWN','stable') )

)

table(deg$g)

up_GSE65194=rownames(deg[deg$g=='UP',])

down_GSE65194=rownames(deg[deg$g=='DOWN',])

save(up_GSE65194,up_GSE45827,up_GSE38959,

down_GSE65194,down_GSE45827,down_GSE38959,

file = 'all_up_down.rdata')

load(file = 'all_up_down.rdata')

library(VennDiagram)

venn.diagram(list( up_GSE65194=up_GSE65194 ,

up_GSE45827=up_GSE45827,

up_GSE38959=up_GSE38959 ),

fill=c("red","green","blue"),

filename="up_VennDiagram.tiff")

venn.diagram(list( down_GSE65194=down_GSE65194 ,

down_GSE45827=down_GSE45827,

down_GSE38959=down_GSE38959 ),

fill=c("red","green","blue"),

filename="down_VennDiagram.tiff")

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

image-20191121212838253

最后友情宣传生信技能树

9人点赞

数据挖掘学习记录

"小礼物走一走,来简书关注我"

还没有人赞赏,支持一下

小梦游仙境

总资产46 (约4.28元)共写了21.3W字获得593个赞共667个粉丝

点击数:0