R语言在基因芯片数据处理中的应用要点

1.R语言安装：官方网站/安装软件。

2. 所需要的软件包：

2.1 affy数据处理相关的程序包

在R中复制source("/biocLite.R")

biocLite("affy")

2.2 热度图相关程序包

Gplots（）：install.packages("gplots")

3.获取基因表达数据

3.1 读取基因芯片数据（cel.files）

the.filter <- matrix(c("CEL file (*.cel)", "*.cel", "All (*.*)", "*.*"), ncol = 2, byrow = T)

cel.files <- choose.files(caption = "Select CEL files", multi = TRUE, filters = the.filter, index = 1)

raw.data <- ReadAffy(filenames = cel.files)

3.2 sampleNames(raw.data)ang #先看看原样品名称的规律

3.3 pat <- ".*MT-([0-9A-Z]+).*" #样品名称查找的正则表达式

sampleNames(raw.data) <- gsub(pat, "\\1", cel.files) #gsub为正则表达式查找函数(samples <- sampleNames(raw.data))

pData(raw.data)$treatment <- rep(c("0h", "1h", "24h", "7d"), each = 2)#确定样品重复数使用rma方法进行预处理：eset.rma <- rma(raw.data)

## Background correcting

## Normalizing

## Calculating Expression

4.计算基因表达量

emat.rma.log2 <- exprs(eset.rma)

class(emat.rma.log2)

head(emat.rma.log2, 1)

# 计算平均值，并做对数转换

results.rma <- data.frame((emat.rma.log2[, c(1, 3, 5, 7)] + emat.rma.log2[, c(2, 4, 6, 8)])/2)

# 计算表达量差异倍数

results.rma$fc.1h <- results.rma[, 2] - results.rma[, 1]

results.rma$fc.24h <- results.rma[, 3] - results.rma[, 1]

results.rma$fc.7d <- results.rma[, 4] - results.rma[, 1]

head(results.rma, 2)

5.T检验

p.value = apply(emat.rma.log2,1, function(x)(t.test(x[7:9], x[10:12])$p.value))

6.导出数据

write.csv(results.rma,file="C:/users/suntao/desktop/data.csv")

7. 选取目的基因

在/index.php上确定探针，选取数据；汇总到excel表格中，保存为csv格式。

8.热度图

cipk=read.csv("c:/users/suntao/desktop/TaCIPK affx arry log.csv")

s(cipk)=cipk$genename

cipk <- cipk[,-1]

cipk_matrix=data.matrix(cipk)

library(gplots)

heatmap.2(cipk_matrix,Rowv=FALSE,Colv=FALSE,col=greenred(75),key=TRUE,keysize=0.8,trace="n one",="none",symkey=FALSE,revC=FALSE,margins=c(10,10),denscol=tracecol,distfun=dist, hclustfun=hclust,dendrogram="none",symm=FALSE)

heatmap.2颜色选择函数col=colorRampPalette(c("black","red"))

用R和BioConductor进行基因芯片数据分析(二)：缺失值填充

以下分析用到的数据可以在这里(/u/308058/blog/raw_data_3_replicates.txt )下载，这个数据来自关于基因对蝴蝶迁移性的研究，样本是20个蝴蝶个体，其中10个是当地固有个体(old)，另外10个是新迁入的个体(new)，old和new个体两两随机配对，分别用不同颜色染料（波长分别为555和647nm）标记后，在同一张基因芯片上杂交；此外，每个基因在每张芯片上都重复点样3次，因此此数据是有3个replicates 及10张芯片的双通道芯片。数据是样点的信号强度值，没有经过标准化处理的。

拿到数据你会看到许多”NA”，这是因为我把缺失的空白值替换成NA了，以便用R进行缺失值填充。

说到缺失值填充，通常有3种方法:

A. 用此基因的平均表达值填充；如果有多张重复芯片，可以取不同芯片上的平均值；对于时间序列芯片，可以通过插值法填充。此方法很简单，也比较常用，但是效果不及下面2种方法

B. 基于SVD（即单值分解)方法的填充：简单地讲，此方法是通过描述基因表达谱的几个基本模式来对缺失值进行填充。

C. 基于KNN（最近邻）方法的填充：此方法是寻找和有缺失值的基因的表达谱相似的其他基因，通过这些基因的表达值（依照表达谱相似性加权）来填充缺失值。KNN法是这3种方法里效果最好的，因此对本数据的缺失值用KNN法填充。

对以上3种方法的比较，这篇paper提供了清晰的说明: Troyanskaya, O., Cantor, M., Sherlock, G., Brown, P., Hastie, T., Tibshirani, R., Botstein, D., and Altman, R. B. (2001), Missing value estimation methods for DNA microarrays, Bioinformatics 17(6):520-525. 其中关于KNN的介绍如下：

The KNN-based method selects genes with expression profiles similar to the gene of interest to impute missing values. If we consider gene A that has one missing value in experiment 1, this method would find K other genes, which have a value present in experiment 1, with expression most similar to A in experiments 2–N (where N is the total number of experiments). A weighted average of values in experiment 1 from the K closest genes is then used as an estimate for the missing value in gene A.

In the weighted average, the contribution of each gene is weighted by similarity of its expression to that of gene A.