Step2-Differential-Expression-Genes

上一篇中做了：GEO数据下载和表达矩阵提取及质控。接下来是差异基因的获得。

一、差异分析

1、表达矩阵

#1、表达矩阵
load('GSE63067_exprSet_eSet.Rdata')
exprSet <- exprSet

#查看内参表达量

exprSet['GAPDH',]

## GSM1539877 GSM1539878 GSM1539879 GSM1539880 GSM1539881 GSM1539882 
##   12.05812   12.19904   11.85656   11.83436   11.82101   12.04682 
## GSM1539883 GSM1539884 GSM1539885 GSM1539886 GSM1539887 GSM1539888 
##   11.90930   12.14537   12.00037   11.54590   12.28213   12.11924 
## GSM1539889 GSM1539890 GSM1539891 GSM1539892 GSM1539893 GSM1539894 
##   12.15863   11.62218   11.67380   12.05812   12.19178   11.78575

boxplot(exprSet, las = 2)

exprSet['ACTB',]

## GSM1539877 GSM1539878 GSM1539879 GSM1539880 GSM1539881 GSM1539882 
##   12.53594   12.26915   11.90930   12.18436   11.91843   12.06563 
## GSM1539883 GSM1539884 GSM1539885 GSM1539886 GSM1539887 GSM1539888 
##   12.05812   11.87208   11.66478   11.95543   11.85656   12.28213 
## GSM1539889 GSM1539890 GSM1539891 GSM1539892 GSM1539893 GSM1539894 
##   12.99842   11.71382   12.08715   11.82830   12.33032   11.57166

boxplot(exprSet, las = 2)

image.png

2、分组矩阵

group <- c(rep("Steatosis",2),rep("Non_alcoholic_steatohepatitis",9),rep("Healthy",7)) 
head(group)

## [1] "Steatosis"                     "Steatosis"                    
## [3] "Non_alcoholic_steatohepatitis" "Non_alcoholic_steatohepatitis"
## [5] "Non_alcoholic_steatohepatitis" "Non_alcoholic_steatohepatitis"

group <- factor(group)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
design

##    Healthy Non_alcoholic_steatohepatitis Steatosis
## 1        0                             0         1
## 2        0                             0         1
## 3        0                             1         0
## 4        0                             1         0
## 5        0                             1         0
## 6        0                             1         0
## 7        0                             1         0
## 8        0                             1         0
## 9        0                             1         0
## 10       0                             1         0
## 11       0                             1         0
## 12       1                             0         0
## 13       1                             0         0
## 14       1                             0         0
## 15       1                             0         0
## 16       1                             0         0
## 17       1                             0         0
## 18       1                             0         0
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"

3、差异比较矩阵

library(limma)

contrast.matrix <- makeContrasts(Non_alcoholic_steatohepatitis - Healthy, 
                                 Steatosis - Healthy, 
                                 Non_alcoholic_steatohepatitis - Steatosis, 
                                 levels=design)
contrast.matrix

##                                Contrasts
## Levels                          Non_alcoholic_steatohepatitis - Healthy
##   Healthy                                                            -1
##   Non_alcoholic_steatohepatitis                                       1
##   Steatosis                                                           0
##                                Contrasts
## Levels                          Steatosis - Healthy
##   Healthy                                        -1
##   Non_alcoholic_steatohepatitis                   0
##   Steatosis                                       1
##                                Contrasts
## Levels                          Non_alcoholic_steatohepatitis - Steatosis
##   Healthy                                                               0
##   Non_alcoholic_steatohepatitis                                         1
##   Steatosis                                                            -1

fit <- lmFit(exprSet, design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2)

allDiff1=topTable(fit2,adjust='fdr',coef=1,number=Inf) 
allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf) 
allDiff3=topTable(fit2,adjust='fdr',coef=3,number=Inf)

write.csv(allDiff1, file = "allDiff1.csv" );
write.csv(allDiff2, file = "allDiff2.csv" );
write.csv(allDiff3, file = "allDiff3.csv" );