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( = cel.files)
3.2 sampleNames(raw.data)ang #先看看原样品名称的规律
7. 选取目的基因
在上确定探针，选取数据；汇总到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"))
中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 expe riment 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.
下面分析数据
首先要安装R（)
然后下载安装叫做impute的package
> source("")
> biocLite("impute")
impute是专门用KNN法进行缺失值填充的R package:
设置好当前工作目录( Windows是在R的菜单栏->文件->改变工作目录…设置，Linux下用setwd()函数）然后在R控制台输入以下代码:
library(impute)
#导入impute package
raw<-read.table('raw_data_3_replicates.txt',header=TRUE)
rawexpr<-raw[,-1]
#移除第一列ID列
if(exists(".Random.seed")) rm(.Random.seed)
#必须，如果没有这句话会出错，原因不知-,-请高手指教
imputed<-impute.knn(as.matrix(rawexpr) ,k = 10, rowmax = 0.5, colmax = 0.8, maxp = 1500, rng.seed=362436069)
#impute.knn() 使用一个矩阵作为第一个参数，其他参数这里使用的是默认值
write.table(imputed$data,file='imputed_data.txt')
#write.table() 把数据保存在当前工作目录下的文件中，文件名用file=’ ‘指定，这一步不是必须的imputeddata<-imputed$data
#imputed$data是在R中储存imputed后的数据的矩阵
现在在R里输入imputed，即填充好的数据矩阵，是不是NA值全都没了？
关于impute package的详细Documentation在
全部数据文件：
from :
用R和BioConductor进行基因芯片数据分析(三)：计算median
接前一篇：
我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。

这一步很简单，就是取这3个值的中位数，即median。

方法很多，在excel中可以用median函数;
在R中以下代码进行操作：
get_median<-function(i,j){
num_vec<-c(imputeddata[i*3-2,j],imputeddata[i*3-1,j],imputeddata[i*3,j])
median(num_vec)
}
#A simple function to calculate median value of three replicates
dimrow<-(dim(imputeddata)[1])/3
mediandata<-matrix(data = NA, nrow =dimrow, ncol = dim(imputeddata)[2], byrow = TRUE, dimnames = NULL)
#Create a blank matrix to store median values
for (i in 1:dimrow){
for (j in 1:dim(imputeddata)[2]){
mediandata[i,j]<-get_median(i,j)
}
}
#Assign median value using the function get_median()
现在我们得到了中位数的数据，储存在mediandata对象里，行数是缺失值填充数据imputeddata的1/3，double check一下：
> dim(imputeddata)
[1] 11571 20
> dim(mediandata)
[1] 3857 20
from:
用R和BioConductor进行基因芯片数据分析(三)：计算median
接前一篇：
我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。

这一步很简单，就是取这3个值的中位数，即median。

方法很多，在excel中可以用median函数;
在R中以下代码进行操作：
get_median<-function(i,j){
num_vec<-c(imputeddata[i*3-2,j],imputeddata[i*3-1,j],imputeddata[i*3,j])
median(num_vec)
}
#A simple function to calculate median value of three replicates
dimrow<-(dim(imputeddata)[1])/3
mediandata<-matrix(data = NA, nrow =dimrow, ncol = dim(imputeddata)[2], byrow = TRUE, dimnames = NULL)
#Create a blank matrix to store median values
for (i in 1:dimrow){
for (j in 1:dim(imputeddata)[2]){
mediandata[i,j]<-get_median(i,j)
}
}
#Assign median value using the function get_median()
现在我们得到了中位数的数据，储存在mediandata对象里，行数是缺失值填充数据imputeddata的1/3，double check一下：
> dim(imputeddata)
[1] 11571 20
> dim(mediandata)
[1] 3857 20
from:
用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化
接前一篇：用R和BioConductor进行基因芯片数据分析(四)：芯片内归一化
上次进行了芯片内的归一化，但是我们的数据来自于10张芯片，为了让这10张芯片之间有可比性，需要进行芯片间归一化。

具体原理就不介绍了。

这里用到Bioconductor的一个package，叫做limma，以及其中的函数normalizeBetweenArrays()
由于normalizeBetweenArrays()需要log intensity或log ratio作为输入，于是先进行log转化：
#log transformation
norm_log<-matrix(data = NA, nrow =dim(normed)[1], ncol = dim(normed)[2], byrow = TRUE, dimnames = NULL)
for (i in 1:dim(normed)[1]){
for (j in 1:dim(normed)[2]){
norm_log[i,j]<-log(normed[i,j])/log(2)
}
}
然后利用函数进行芯片间归一化：
source("")
biocLite("limma")
library(limma)
norm_log_btw_array<-normalizeBetweenArrays(norm_log,method='scale') normalizeBetweenArrays()函数有许多方法，具体请看帮助。

下面看看效果吧
plot(density(norm_log[,1]),ylim=c(0,1.35),xlab='log intensity')
for (i in 2:20){
lines(density(norm_log[,i]),type='l')
}
lines(density(norm_log_btw_array[,1]),type='l',col='green')
for (i in 2:20){
lines(density(norm_log_btw_array[,i]),type='l',col='green')
}
text(1.5,c(0.8,1.0),labels=c('BEFORE normalization','AFTER normalization'),col=c('black','green'))
用R和BioConductor进行基因芯片数据分析(六)：差异表达基因
接前一篇：用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化
经过一系列的预处理，包括缺失值填充，中位数计算以及归一化，我们的数据终于可以用啦。

下面我们就来分析一下new population和old population的个体是否有差异表达基因。

判断一个基因是否差异表达有许多方法，最早使用的就是看log ratio的绝对值是否大于2,这种方法早已废弃。

下一个想到的也许是t-test，诚然t-test可以统计地判断一个基因是否差异表达，但是对于有数千数万基因的芯片来说，会有很高的错误发现率（False Discovery Rate, FDR)，如果p value < 0.05,则10000个基因里有500个基因实际没有差异表达而被误认为是差异表达。

因此t-test方法需要改进。

于是Westfall & Young (1993) 提出了Step-down maxT and minP multiple testing procedures，大意就是比较几个group间有没有差异基因表达，就通过随机置换这些group的标记，相当于随机互换group的成员，模拟一个空分布(null distribution），以此计算调整后的p value，这个方法可以极大的减小Family-wise Error Rate (FWER).
以下分析就使用Step-down maxT and minP multiple testing procedures，需要求助于Bioconductor的multtest package的mt.maxT()函数。

具体用法可通过help(mt.maxT)查看。

source("")
biocLite("multtest")
library(multtest)
classlabel<-c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)
maxttt<-mt.maxT(norm_log_btw_array,classlabel,B=100000)
默认随机置换次数B=10000，对于microarray来说B应该比10000大很多，所以这里取B=100000
以下是画图：
rawp<-maxttt$rawp[order(maxttt$index)]
plot(sort(rawp),type='p',col=1,ylim=c(-0.05,1.00),ylab='p value')
lines(sort(maxttt$adjp),type='p',col='red')
#min adj-p: sort(maxttt$adjp)[1] 0.0163
#rawp: >sort(rawppp)[170] [1] 0.0493 > sort(rawppp)[171] [1] 0.0502 170个raw p小于0.05 abline(h=0.05,col='blue')
text(1000,c(0.6,0.7),labels=c('raw p-value','adjusted p-value'),col=c('black','red'))
text(1000,0.08,labels='p=0.05',col='blue')
可见调整后只有一个基因的p value小于0.05，而未调整的有170个基因的p value小于0.05，可以说虽然此方法降低了错误发现率，但是也导致了很高的False negative.
此外可以考虑使用multtest package的mt.rawp2adjp()函数，这个函数可以通过”Bonferroni”, “Holm”, “Hochberg”, “SidakSS”, “SidakSD”, “BH”, “BY”等方法调整p value，不过对我们的数据来说都过于严格了。

procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY")
adjps<-mt.rawp2adjp(rawp,procs)
plot(sort(adjps$adjp[,1]),ylab='p value')
for (i in 2:8){
points(sort(adjps$adjp[,i]),col=i)
}
abline(h=0.05,col='blue')
text(1000,0.08,labels='p=0.05',col='blue')
因此可以考虑不这么严格的SAM (Significance Analysis of Microarrays)分析。

有兴趣的请看参考资料。

参考资料：
课堂讲义：Differential expression. Identifying differentially expressed genes — notions on multiple testing and p-value adjustments.
Dudoit, S., Yang, Y.H., Speed, T.P., and Callow, M.J. (2002), Statistical methods for identifying differentially expressed genes in replicated cDNA microarray
experiments, Statistica Sinica 12(1):111-139.
Dudoit S., Shaffer J.P., Boldrick J.C. (2003). Multiple hypothesis testing in microarray experiments, Statistical Science, 18(1): 71-103.
Efron B., Tibshirani, R., Storey J.D., and Tusher V. (2001), Empirical Bayes analysis of a microarray experiment, Journal of the American Statistical Association 96:1151-1160.
SAM (Significance Analysis of Microarrays)相关：Tusher, V.G., Tibshirani, R., and Chu, G. (2001) Significance analysis of microarrays applied to the ionizing radiation
response, PNAS 98:5116-5121.
SAM 网站
FDR相关：Storey J.D. (2002), A direct approach to false discovery rates, JRSS-B 64(3):479-498.
from :
安装R package的2种方法
安装R语言的包的方法：
1. 自动安装（在线安装）
在R的控制台,输入
install.packages(“stepNorm”,contriburl=””, dependencies = TRUE)
若要指定安装目录(e.g. “mydir”)，则输入
install.packages(“stepNorm”,contriburl=””,lib=”mydir”)
这种方法可能找不到需要的package，那么可以用方法2
2. 手动安装（离线安装）
Windows:
下载package.zip文件
打开R的菜单栏->Packages->“Install package from local zip file…”
选择package.zip文件
Linux上安装R包（离线安装）:
下载package.tar.gz文件
在Shell终端（注意不是R）输入:
sudo R CMD INSTALL package.tar.gz
注意：需要sudo权限才能安装。

否则提示：
username is not in the sudoers file. This incident will be reported
如何把sra格式转成fastq格式(fq格式)
sra是NCBI 推出的存储高通量数据的格式，而平常我们工作用得多是fastq格式。

如果需要把sra 转成fastq，从
下载相应的软件。

或者下载最新的source code，在服务器上用make 编译。

然后使用如下命令行：
sra_sdk-2.0.0rc1/linux/rel/gcc/x86_64/bin/fastq-dump -A SRR034580 -D SRR034580.sra 这样就可以很简单的把sra格式转成fastq格式了。

REF:
一、R语言相关资料
1、R语言
2、R语言课件（复旦大学）
3、上传：GLM课件（R语言）
4、利用R语言实现微阵列数据分析-聚类分析
5、Medline 文献挖掘的开放资源库----MedlineR
二、R语言相关问题
1、R 语言,请帮忙
2、请教:有用过R语言和fbioconductor的吗
3、求助：perl程序（加一点统计知识）
4、请教诸位高手关于单倍型分析的一些问题
5、perl 求助一程序。

6、[求助]关于SAM 谁能帮我讲解一下SAM方法？我不是很懂
R语言相关网站(超有用资源!!!!!!!!!!!!!):
R-FAQ的版本在CRAN网站定期更新:
有关R或者统计方法的出版物:
R网站地址清单:
提交问题"公告指南" 电子杂志R News:
R Graphical Manuals:。