beat365官网_beat365亚洲官方网站

  • beat365官网
  • |
  • beat365官网生物
  • |
  • beat365官网诊断
  • |
  • beat365官网检验


beat365
  • 免费服务热线
  • 400-065-6886
  • 电话:86(0)512-6295 9990
  • 传真:86(0)512-6295 9995
新闻中心

如何对GEO数据进行差异分析

发稿时间:2020-06-16来源:beat365官网生物





GEO数据库目前收录了4348个数据集记录,包含人的数据1772个,小鼠的数据1642个,大鼠的数据360个,其中属于组织样品有1183个,细胞品系有857个。下面,小编就跟大家详细讲解如何利用GEO表达谱数据进行差异表达分析,期待您的评论沟通和转发哦~

beat365



GEO数据下载

beat365


打开NCBI官网,选择GEO DataSets,这里我们随便搜一个转录组的数据,如上图所示,由于前面几个GSE所提供的表达量文件不规范,这里我们选择登录号为GSE132287,点击进入。 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132287


beat365


下载文件打开后,如图所示:

beat365



数据准备

1、表达量数据读取

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. 分组数据读取

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数据分析

1、安装和加载DESeq2包

In [5]:

install.packages('BiocManager')
BiocManager::install('DESeq2')

In [6]:

library(DESeq2)


2、DESeq2分析

构建数据集并标准化数据集

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

表达量数据归一化

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倍数变化和校正前后p值等信息。contrast参数用于指定比较的分组顺序,即谁相对于谁的表达量上调/或下调;pAdjustMethod设定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


合并数据

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)


往期相关链接:

如何使用Rstudio练习R基础教程

R相关软件及R包安装

【绘图进阶】之交互式可删减分组和显示样品名的PCA 图(三)

【绘图进阶】之绘制PCA biplot图(二);

【进阶篇绘图】之带P值的箱体图、小提琴图绘制(一);

3分钟学会CHIP-seq类实验测序数据可视化 —IGV的使用手册;

10分钟搞定多样性数据提交,最快半天内获取登录号,史上最全的多样性原始数据提交教程

【零基础学绘图】之绘制venn图(五);

【WGS服务升级】人工智能软件SpliceAI助力解读罕见和未确诊疾病中的非编码突变

【零基础学绘图】之绘制barplot柱状图图(四)

【零基础学绘图】之绘制heatmap图(三)

20分钟搞定GEO上传,史上最简单、最详细的GEO数据上传攻略

【零基础学绘图】之绘制PCA图(二)

【零基础学绘图】之alpha指数箱体图绘制(一)

beat365

如果您对本文案介绍的方法或代码有疑问,

请扫码添加QQ群沟通

【本群将为大家提供】

分享生信分析方案

提供数据素材及分析软件支持

定期开展生信分析线上讲座

beat365

QQ号:1040471849

Copyright ? 2012-2020 eat365官网,beat365亚洲官方网站   All Rights Reserved    苏ICP备17064027号-1
XML 地图 | Sitemap 地图