如何對GEO數(shù)據(jù)進行差異分析

發(fā)稿時間：2020-06-16來源：天昊生物

GEO數(shù)據(jù)庫目前收錄了4348個數(shù)據(jù)集記錄，包含人的數(shù)據(jù)1772個，小鼠的數(shù)據(jù)1642個，大鼠的數(shù)據(jù)360個，其中屬于組織樣品有1183個，細胞品系有857個。下面，小編就跟大家詳細講解如何利用GEO表達譜數(shù)據(jù)進行差異表達分析，期待您的評論溝通和轉(zhuǎn)發(fā)哦~

一GEO數(shù)據(jù)下載

打開NCBI官網(wǎng)，選擇GEO DataSets，這里我們隨便搜一個轉(zhuǎn)錄組的數(shù)據(jù)，如上圖所示，由于前面幾個GSE所提供的表達量文件不規(guī)范，這里我們選擇登錄號為GSE132287，點擊進入。 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132287

下載文件打開后，如圖所示：

二數(shù)據(jù)準備

1、表達量數(shù)據(jù)讀取

In [1]:

gene <- read.table('GSE132287_Gene-count-table.xls',header = T, row.names = 1, sep = 't', check.names = FALSE)

head(gene,10)

Out[1]:

 gene_name gene_type chr start end strand length MDA-A1 MDA-A2 MDA-A3 MDA-G1 MDA-G2 MDA-G3

ENSG00000000003.14_2 TSPAN6 protein_coding chrX 99882106  99894988 - 4535  2704  3694  2946  3204  2435  3498

ENSG00000000005.5_2 TNMD protein_coding chrX 99839799  99854882  +  1610  0  0  0  0  0  0

ENSG00000000419.12_2 DPM1 protein_coding chr20 49551404  49575092 - 1207  4045  5254  4031  4740  3872  4867

ENSG00000000457.13_2 SCYL3 protein_coding chr1 169818772  169863408 - 6883  1352  1849  1431  1546  1291  1752

ENSG00000000460.16_4 C1orf112 protein_coding chr1 169631245  169823221  +  5967  3490  4828  4071  3885  3376  4586

ENSG00000000938.12_2 FGR protein_coding chr1 27938575  27961788 - 3474  1  0  0  1  0  3

ENSG00000000971.15_2 CFH protein_coding chr1 196621008  196716634  +  8145  1620  2346  2001  2148  1863  2507

ENSG00000001036.13_2 FUCA2 protein_coding chr6 143815948  143832827 - 2793  9471  12928  10236  11808  10084  13539

ENSG00000001084.10_3 GCLC protein_coding chr6 53362139  53481768 - 8463  2461  3384  2572  2761  2392  3046

ENSG00000001167.14_2  NFYA  protein_coding  chr6  41040684  41067715  +  3811  4549  4770  3828  4146  4264  4997

2. 分組數(shù)據(jù)讀取

In [2]:

info= read.table('GSE132287_sample_group.txt',header = F,col.names = c('sample','group'))

info

Out[2]:

sample group

MDA-A1  MDA-A

MDA-A2  MDA-A

MDA-A3  MDA-A

MDA-G1  MDA-G

MDA-G2  MDA-G

MDA-G3  MDA-G

In [3]:

df = gene[as.character(info$sample)]

head(df)

Out[3]:

 MDA-A1 MDA-A2 MDA-A3 MDA-G1 MDA-G2 MDA-G3

ENSG00000000003.14_2 2704  3694  2946  3204  2435  3498

ENSG00000000005.5_2 0  0  0  0  0  0

ENSG00000000419.12_2 4045  5254  4031  4740  3872  4867

ENSG00000000457.13_2 1352  1849  1431  1546  1291  1752

ENSG00000000460.16_4 3490  4828  4071  3885  3376  458

6ENSG00000000938.12_2  1  0  0  1  0  3

In [4]:

coldata <- data.frame(group = info$group )

coldata

Out[4]:

group

MDA-A

MDA-A

MDA-A

MDA-G

MDA-G

MDA-G

三Deseq數(shù)據(jù)分析

1、安裝和加載DESeq2包

In [5]:

install.packages('BiocManager')

BiocManager::install('DESeq2')

In [6]:

library(DESeq2)

2、DESeq2分析

構(gòu)建數(shù)據(jù)集并標準化數(shù)據(jù)集

In [7]:

dds <- DESeqDataSetFromMatrix(df, DataFrame(coldata), design= ~ group )

dds <- DESeq(dds,betaPrior=FALSE)

dds

Out[7]:

class: DESeqDataSet

dim: 60461 6

metadata(1): version

assays(4): counts mu H cooks

rownames(60461): ENSG00000000003.14_2 ENSG00000000005.5_2 ...ENSG00000284747.1_1 ENSG00000284748.1_1

rowData names(22): baseMean baseVar ... deviance maxCooks

colnames(6): MDA-A1 MDA-A2 ... MDA-G2 MDA-G3

colData names(2): group sizeFactor

表達量數(shù)據(jù)歸一化

In [8]:

df <- as.data.frame(counts(dds, normalized=TRUE))

In [9]:

df['MDA-A_mean'] = apply(

df[as.character(info[info$group=='MDA-A',1])],1,mean)

df['MDA-G_mean'] = apply(

df[as.character(info[info$group=='MDA-G',1])],1,mean)

df['gene'] = rownames(df)head(df)

Out[9]:

 MDA-A1 MDA-A2 MDA-A3 MDA-G1 MDA-G2 MDA-G3 MDA-A_mean MDA-G_mean gene

ENSG00000000003.14_2 2994.378835 3147.518 3196.921 3101.6740411 2819.888 2992.716062 3112.9393704 2971.425929 

ENSG00000000003.14_2ENSG00000000005.5_2 0.000000 0.000 0.000 0.0000000 0.000 0.000000 0.0000000 0.000000 

ENSG00000000005.5_2ENSG00000000419.12_2 4479.386978 4476.735 4374.335 4588.6188998 4484.027 4163.964858 4443.4855603 4412.203499 

ENSG00000000419.12_2ENSG00000000457.13_2 1497.189418 1575.463 1552.883 1496.6254893 1495.062 1498.924683 1541.8452966 1496.870591 

ENSG00000000457.13_2ENSG00000000460.16_4 3864.786292 4113.757 4417.742 3760.9249843 3909.627 3923.555134 4132.0948079 3864.702246 

ENSG00000000460.16_4ENSG00000000938.12_2  1.107389  0.000  0.000  0.9680631  0.000  2.566652  0.3691295  1.178238  ENSG00000000938.12_2

差異分析

通過result()可獲得最終計算的log2倍數(shù)變化和校正前后p值等信息。contrast參數(shù)用于指定比較的分組順序，即誰相對于誰的表達量上調(diào)/或下調(diào)；pAdjustMethod設(shè)定p值校正方法；alpha為顯著性水平，這里0.05為校正后p值小于0.05即為顯著。In [10]:

res <- results(dds, contrast = c('group', 'MDA-A', 'MDA-G'), pAdjustMethod = 'fdr', alpha = 0.05)

res = as.data.frame(res)

head(res)

Out[10]:

baseMean log2FoldChange lfcSE stat pvalue padj

ENSG00000000003.14_2 3042.1826497 0.06681719 0.07222184 0.9251661 0.3548795 0.7074537

ENSG00000000005.5_2 0.0000000 NA NA NA NA NAE

NSG00000000419.12_2 4427.8445298 0.01055749 0.06762388 0.1561207 0.8759379 0.9691177

ENSG00000000457.13_2 1519.3579439 0.04308650 0.07852322 0.5487103 0.5832043 0.8550993

ENSG00000000460.16_4 3998.3985272 0.09654288 0.07227280 1.3358121 0.1816107 0.5297536

ENSG00000000938.12_2  0.7736839  -1.73212025  2.76415576  -0.6266363  0.5308977  NA

In [11]:

res['type']='Not DEG'

res[which(res$log2FoldChange >= 1 & res$pvalue < 0.05),'type'] <- 'Up'

res[which(res$log2FoldChange <= 1 & res$pvalue < 0.05),'type'] <- 'Down'

res['gene'] = rownames(res)head(res)

Out[11]:

 baseMean log2FoldChange lfcSE stat pvalue padj type gene

ENSG00000000003.14_2 3042.1826497 0.06681719 0.07222184 0.9251661 0.3548795 0.7074537 Not DEG E

NSG00000000003.14_2ENSG00000000005.5_2 0.0000000 NA NA NA NA NA Not DEG ENSG00000000005.5_2

ENSG00000000419.12_2 4427.8445298 0.01055749 0.06762388 0.1561207 0.8759379 0.9691177 Not DEG 

ENSG00000000419.12_2ENSG00000000457.13_2 1519.3579439 0.04308650 0.07852322 0.5487103 0.5832043 0.8550993 Not DEG

ENSG00000000457.13_2ENSG00000000460.16_4 3998.3985272 0.09654288 0.07227280 1.3358121 0.1816107 0.5297536 Not DEG 

ENSG00000000460.16_4ENSG00000000938.12_2  0.7736839  -1.73212025  2.76415576  -0.6266363  0.5308977  NA  Not DEG 

ENSG00000000938.12_2

合并數(shù)據(jù)

In [12]:

result_merge = merge(df,res,by = 'gene')

result_merge = result_merge[order(result_merge$pvalue),]

head(result_merge)

Out[12]:

gene MDA-A1 MDA-A2 MDA-A3 MDA-G1 MDA-G2 MDA-G3 MDA-A_mean MDA-G_mean baseMean log2FoldChange lfcSE stat pvalue padj type1875

 ENSG00000090339.8_2 4852.57694 4787.73773 5220.7701 27261.624 26204.689 25266.976 4953.69492 26244.430 15599.062 -2.405646 0.06409772 -37.53091 0.000000e+00  0.000000e+00 Down11138 

ENSG00000163739.4_2 193.79301 177.22895 210.5237 3221.714 3519.359 3184.359 193.84854 3308.477 1751.163 -4.097951 0.10170786 -40.29139 0.000000e+00  0.000000e+00 Down11703 

ENSG00000165795.23_3 17.71822 28.11805 41.2366 45895.870 41106.667 43346.472 29.02429 43449.669 21739.347 -10.548113 0.16963433 -62.18148 0.000000e+00  0.000000e+00 Down12624 

ENSG00000169429.10_2 2099.60883 2407.92789 2647.8235 28439.757 28773.277 27150.899 2385.12008 28121.311 15253.215 -3.559289 0.07909170 -45.00205 0.000000e+00  0.000000e+00 Down10535

ENSG00000160710.15_3 85045.23143 85230.08188 83116.6997 21604.263 22847.460 22927.901 84464.00435 22459.875 53461.939 1.911000 0.05663594 33.74183 1.409012e-249 4.237464e-246 Up11364  

ENSG00000164400.5_2  759.66860  834.16893  959.2934  6727.070  5916.553  6392.674  851.04366  6345.432  3598.238  -2.898827  0.08922127  -32.49032  1.460968e-231  3.661428e-228  Down

寫入文件保存

In [13]:

write.csv(result_merge,file = 'Differential_Expression_Genes_Summary.csv', quote=F,row.names =F)

往期相關(guān)鏈接：

如何使用Rstudio練習(xí)R基礎(chǔ)教程；

R相關(guān)軟件及R包安裝；

【繪圖進階】之交互式可刪減分組和顯示樣品名的PCA 圖（三）；

【繪圖進階】之繪制PCA biplot圖(二)；

【進階篇繪圖】之帶P值的箱體圖、小提琴圖繪制（一）；

3分鐘學(xué)會CHIP-seq類實驗測序數(shù)據(jù)可視化 —IGV的使用手冊；

10分鐘搞定多樣性數(shù)據(jù)提交，最快半天內(nèi)獲取登錄號，史上最全的多樣性原始數(shù)據(jù)提交教程；

【零基礎(chǔ)學(xué)繪圖】之繪制venn圖（五）；

【W(wǎng)GS服務(wù)升級】人工智能軟件SpliceAI助力解讀罕見和未確診疾病中的非編碼突變；

【零基礎(chǔ)學(xué)繪圖】之繪制barplot柱狀圖圖（四）；

【零基礎(chǔ)學(xué)繪圖】之繪制heatmap圖（三）；

20分鐘搞定GEO上傳，史上最簡單、最詳細的GEO數(shù)據(jù)上傳攻略；

【零基礎(chǔ)學(xué)繪圖】之繪制PCA圖（二）；

【零基礎(chǔ)學(xué)繪圖】之a(chǎn)lpha指數(shù)箱體圖繪制（一）；

如果您對本文案介紹的方法或代碼有疑問，

請掃碼添加QQ群溝通

【本群將為大家提供】

分享生信分析方案

提供數(shù)據(jù)素材及分析軟件支持

定期開展生信分析線上講座

QQ號：1040471849