affy 包处理affymetrix 表达谱芯片

在正式开始操作前做了一些资料的收集，正好搜到了18级学姐的CSDN文章，从她的作业流程中参考了很多↓

用Bionconductor的affy包处理.cel文件_affy包怎么用

还记得大二学统计的时候，段哥苦口婆心劝我们多看help，这次作业Help还help了我挺多的

整个作业的流程我用到的三个包

library(affy)
library(GEOquery)
library(dplyr)

1. 从Gene Expression Omnibus (GEO) 查询GSE65496，并下载GSE65496_family.soft和原始数据(.cel格式)

老老实实查询的话

红框是要求的文件的信息，蓝色是重要信息

当然也可以用GEOquery包里的函数直接抓取的
gteGEPSuppFiles()
区别呢就是到后面read那一步啊前者我把最外层的tar解了之后放在workspace里面，给个路径就行（我自己的做法），后者则需要在获取tar里信息的时候再给个路径

2. 仔细研读soft格式文件

在GEO官网查一查SOFT文件是个什么格式

纯文本，\t是分隔符

其实也可以直接用excel或者wps表格打开

打开慌神了，啥玩意(是一些其他信息)

好的，到这里才是这次实验需要用到的信息

3. 试着用Bioconductor的affy包下函数处理.cel文件 (可以网页搜索，或者上Bioconductor官网查询)，得到探针表达值

untar("GSE65496_RAW.tar")
raw <- ReadAffy(celfile.path = "GSE65496_RAW")

关于ReadAffy()这个函数，看affy包内容的时候其实没看到，而是read.affybatch()。

看得出来我们就是要一个读入cel文件的函数

在help文档里，我们可以得知ReadAffy()已经把read.affybatch()warp了。

ReadAffy是read.affybatch的一个包装器，允许用户使用部件读取phenoData、MIAME信息和CEL文件。人们还可以定义文件来读取PhenoData和MIAME信息。

不用把套娃的cel.tar一个一个解开了，一步到位

4. 下载该数据对应的平台数据，通过文本处理，将探针和gene symbol对应，得到gene symbol的表达值

质控

①灰度图

很简单的一步，但是耗时挺长。如果是par()成4行2列，可能会报figure margin too large的错，网上查找了一下相关的解决方法 ↓

R绘图： figure margins too large错误_figure margins

jpeg("GSE65496_grayscale.jpg", width=200*2, height=200*4)
par(mfcol=c(4,2))
for(i in 1:8){
  image(raw[,i])
}
dev.off()

GSE65496_grayscale.jpg

②RNA降解信息

阅读包说明文档

也是一个自带函数

但是是不提供注释信息的，所以又加了个框框

par(mfcol=c(1,1))
RNAdeg <- AffyRNAdeg(raw)
plotAffyRNAdeg(RNAdeg,col=topo.colors(8),lwd=2)
legend("topleft",sampleNames(raw),col=topo.colors(8),lwd=2,inset = 0.05,cex=0.3)

GSE65496_degregation.jpg

结果

背景矫正和归一化(mas5或rma方法)

两种方法就是函数名 mas5() rma()
原理上的区别的话↓

MAS5算法是由Affymetrix公司设计开发的，其对于芯片背景噪音的校正采用了比较直接的策略，即通过计算PM[i]-MM[ii]或PM/MM实现信噪分离。理论上讲，MM探针只有非特异杂交信号，不会有特异性杂交，因此MM探针的信号值应始终小于其对应的PM探针信号值。但是，在实际分析过程中，经常发现甚至多达三分之一的MM探针信号值高于对应PM探针信号值。
RMA算法并不直接利用MM的信号去除背景噪音，而是基于20组探针的信号分布采用随机模型来评估表达值，对于低噪号的实验有更好的适用性。

我用下来RMA算法是明显的更快速的，不过两个结果都可以放上来看一下(两种函数的结果都是list，不可以直接write.table()，我试过了。用一下expr()再输出)

msa5

可以看到我苦苦探索的过程😂

#write.csv(data.mas5,"GSE65496_mas5.csv")
#csv长度不够用
#write.table(data.mas5,"GSE65496_mas5.txt",sep='\t')
#输出结果不尽如人意，看了看老师的输出结果，可能之后要稍微修改一下
#好吧 本来像个list 有专门的提取的
expr.mas5 <- exprs(data.mas5)
write.table(expr.mas5,"GSE65496_mas5.txt",sep='\t')

mas5结果：综合考虑pm和mm信号，expression没有取对数

rma

data.rma <- rma(raw)
#write.table(data.rma,"GSE65496_rma.txt",sep='\t')
expr.rma <- exprs(data.rma)
write.table(expr.rma,"GSE65496_rma.txt",sep='\t')

rma结果：RMA只使用pm信号，expression已经对数处理过

探针号和gene symbol对应，得到gene表达值

如果前面硬刚不用GEOqueory包还可以，这一步还是要用啊(这里的GPL570就是前面下载的SOFT文件，看我一开始geo页面的蓝框，看到那个platform了没有？)

GPL570 <- getGEO("GPL570")
names <- Table(GPL570)[c(1,11)]

这里table()函数是不好使的，但是table输了一半的时候看到有个Table()函数，来自GEOqueory。

剩下就是一些非常简单但是绕来绕去的提取、连接工作了

mas5

expr.mas5 <-as.data.frame(expr.mas5)
expr.mas5$ID <- rownames(expr.mas5)
processed.mas5 <- inner_join(names,expr.mas5,by="ID")
processed.mas5 <- processed.mas5[,-1]
write.table(processed.mas5,"GSE65496_processed_mas5.txt",sep="\t")

和前面没啥区别，就是探针名字变基因名了

rma

expr.rma <-as.data.frame(expr.rma)
expr.rma$ID <- rownames(expr.rma)
processed.rma <- inner_join(names,expr.rma,by="ID")
processed.rma <- processed.rma[,-1]
write.table(processed.rma,"GSE65496_processed_rma.txt",sep="\t")