基因表达分析

上善若水

求大家来帮我解答下  R语言这样做的对吗

1. library(hgu133a.db);
   
2. 
   
3. oriInData<-read.table("D:/实验1/GSE15852_series_matrix.txt",header=TRUE,sep="\t")
   
4. for(i in 2:length(oriInData[1,])){
   
5. oriInData[,i]<-log2(oriInData[,i])
   
6. }
   
7. label<-read.table("D:/实验1/classlable15852.txt")
   
8. ##############################################################
   
9. ###############处理数据###########################3
   
10. ###############################################################
    
11. x<-hgu133aENTREZID;
    
12. probeGene<-toTable(x);
    
13. #Normally, len should be equal to the number of rows of relationship inforinDataion.
    
14. locations<-match(as.character(oriInData[,1]),probeGene$probe_id);
    
15. probegene<-probeGene$gene_id[locations];
    
16. oriInData[,length(oriInData[1,])+1]<-probegene
    
17. oriInData<-na.omit(oriInData)
    
18. gene<-oriInData[,length(oriInData[1,])];
    
19. probe<-unique(gene)
    
20. inData<-matrix(0,nrow=length(probe), ncol=length(oriInData[1,])-2)
    
21. for(i in 1:length(probe)){
    
22.     if(length(which(gene==probe))==1){
    
23.      inData[i,]<-as.numeric(oriInData[which(gene==probe),2:(length(oriInData)-1)])
    
24.      }
    
25.      else if(length(which(gene==probe))>1){
    
26.      inData[i,]<-as.numeric(colMeans(oriInData[which(gene==probe),2:(length(oriInData)-1)]))
    
27.      }
    
28. }     
    
29. inData<-as.data.frame(inData,row.names=probe)
    
30. names(inData)<-names(oriInData)[2:(length(oriInData[1,])-1)]
    
31. ###############处理类标签######################################
    
32. b<-as.vector(label[1,])
    
33. b[which(b=="normal breast tissue")]<-rep(0,length(which(b=="normal breast tissue")))
    
34. b[which(b=="1breast tumor tissue")]<-rep(1,length(which(b=="breast tumor tissue")))
    
35. classLabel<-as.numeric(b)
    
36. #########################计算T-test函数######################################     
    
37. calT <- function(inData, classLabel) {
    
38.         m <- length(ind1 <- which(classLabel == 1))
    
39.         n <- length(ind2 <- which(classLabel == 0))
    
40.         inData1 <- inData[, ind1]
    
41.         inData2 <- inData[, ind2]
    
42.         rmean1 <- rowMeans(inData1)
    
43.         rmean2 <- rowMeans(inData2)
    
44.         ss1 <- rowSums((inData1 - rmean1)^2)
    
45.         ss2 <- rowSums((inData2 - rmean2)^2)
    
46.         tt <- (m + n - 2)^0.5 * (rmean2 - rmean1)/((1/m + 1/n) *
    
47.             (ss1 + ss2))^0.5
    
48.         return(list(T = tt, df = m + n - 2))
    
49.     }
    
50. deGene<-calT(inData,classLabel)[[1]];
    
51. ####################################################################

复制代码