HTA2.0芯片数据分析比较麻烦,我也尝试给你一些思路

  • 2019 年 10 月 7 日
  • 笔记

表达芯片数据处理教程,早在2016年我就系统性整理了发布在生信菜鸟团博客:http://www.bio-info-trainee.com/2087.html 配套教学视频在B站:https://www.bilibili.com/video/av26731585/ 代码都在:https://github.com/jmzeng1314/GEO

  • 但并不是万能的,即使是表达芯片,也有一些是我的知识盲区,主要是因为没有时间去继续探索整理写教程,毕竟大家都不支持我写教程,没有人打赏或者留言鼓励我!

不过最近学徒问到了[HTA-2_0] Affymetrix Human Transcriptome Array 2.0芯片的事情,其实挺麻烦的,首先需要搞清楚下面3个平台的差异:

  • GPL17586 [HTA-2_0] Affymetrix Human Transcriptome Array 2.0 [transcript (gene) version]
  • GPL19251 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array [probe set (exon) version]
  • GPL16686 [HuGene-2_0-st] Affymetrix Human Gene 2.0 ST Array [transcript (gene) version]

如果是芯片cel原始数据处理

那么看我的教程 你要挖的公共数据集作者上传了错误的表达矩阵肿么办(如何让高手心甘情愿的帮你呢?) 里面有提到:

# BiocManager::install(c( 'oligo' ),ask = F,update = F)  library(oligo)  # BiocManager::install(c( 'pd.hg.u133.plus.2' ),ask = F,update = F)  library(pd.hg.u133.plus.2)    dir='~/Downloads/GSE84571_RAW/'    od=getwd()    setwd(dir)    celFiles <- list.celfiles(listGzipped = T)    celFiles    affyRaw <- read.celfiles( celFiles )    setwd(od)    eset <- rma(affyRaw)    eset    # http://math.usu.edu/jrstevens/stat5570/1.4.Preprocess_4up.pdf    save(eset,celFiles,file = f)    # write.exprs(eset,file="data.txt")  

当然了,你没有R基础是看不懂的哈。自行去B站补充下面的视频知识点:

如果是芯片注释到基因

也是很简单

  library(GEOquery)    #Download GPL file, put it in the current directory, and load it:    gpl <- getGEO('GPL19251', 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)]    head(probe2gene)    library(stringr)    probe2gene$symbol=trimws(str_split(probe2gene$gene_assignment,'//',simplify = T)[,2])    plot(table(table(probe2gene$symbol)),xlim=c(1,50))    head(probe2gene)    save(probe2gene,file='probe2gene.Rdata')  

不同平台, 就替换GPL,然后自行看输入输出判断一下即可,也需要R基础啦。

如果是差异分析

当然了也可以直接修改我的代码,不过这个HTA2.0芯片比较麻烦的在于,基于exon和基于transcript的有一点点区别,因为基于exon理论上可以获取不同剪切体的表达量差异的,虽然实际上很少有人愿意去探索。

发表在 BMC Genomics. 2017; 文章 RNA sequencing and transcriptome arrays analyses show opposing results for alternative splicing in patient derived samples 就提到过。

更麻烦的是,因为这个芯片记录35万个外显子,这样矩阵就非常大, 都可以根据illumina的450K甲基化芯片媲美啦!

可以看到,大量的基因有着1-30个探针,如下:

所以有文章采取非常特殊的差异分析策略:

很大概率上,你是看不懂的

最后留一个作业

大家可以拿 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE118222 数据集,测试看看,走我给大家的代码,得到PCA,热图,火山图,看看是否合理。

GSM3321243    LNCaP EtOH rep 1  GSM3321244    LNCaP EtOH rep 2  GSM3321245    LNCaP DHT rep 1  GSM3321246    LNCaP DHT rep 2  GSM3321247    LNCaP ABT-DHT rep 1  GSM3321248    LNCaP ABT-DHT rep 2  GSM3321249    C4-2 DMSO rep 1  GSM3321250    C4-2 DMSO rep 2  GSM3321251    C4-2 ABT rep 1  GSM3321252    C4-2 ABT rep 2