|
歡迎來到醫(yī)科研，這里是白介素2的讀書筆記，跟我一起聊臨床與科研的故事, 生物醫(yī)學(xué)數(shù)據(jù)挖掘，R語言，TCGA、GEO數(shù)據(jù)挖掘。
在數(shù)據(jù)分析過程中，尤其是在做基因篩選時，常會應(yīng)用到批量篩選，這也是應(yīng)用R語言分析數(shù)據(jù)的優(yōu)勢之一。在這一點(diǎn)上往往在線工具不能提供這樣的功能。
我們來構(gòu)建一個虛擬數(shù)據(jù)，來完成基因之間的批量分析，并導(dǎo)出結(jié)果 先產(chǎn)生一個矩陣 data<-matrix(rnorm(120),nrow=20,ncol = 6)
模仿一個矩陣數(shù)據(jù) rownames(data)<-paste("gene",1:20,sep = "")
colnames(data)<-paste("sample",1:6,sep ="" )
head(data)
## sample1 sample2 sample3 sample4 sample5 sample6
## gene1 0.6713704 0.9258705 0.6584580 0.87819034 1.3191532 -0.18602102
## gene2 -0.3275014 1.2083650 0.1699489 0.54841864 0.3195813 -0.04661051
## gene3 1.3187930 0.8165928 -0.6706707 0.11174014 -0.1358276 -0.43858189
## gene4 -0.8040756 0.3348103 -1.5573817 -0.45613542 -0.3221384 -0.59782362
## gene5 -0.9396341 0.7283671 -1.3389603 -0.50145046 -0.2183778 0.88841818
## gene6 1.1326463 -1.5588096 0.3842358 0.06078977 -1.0196727 0.61900583
實(shí)現(xiàn)批量完成基因間的相關(guān)系數(shù)計算 首先明確實(shí)現(xiàn)這個目的可以使用for循環(huán) ### 構(gòu)建好需要得出的結(jié)果表，包括基因名，相關(guān)系數(shù)，Pvalue
在使用for循環(huán)前，需要先來思考for循環(huán)的三個部分
首先創(chuàng)建空向量 gene_name1<-c()##也可用vector
gene_name2<-c()
cor_r<-c()
pvalue<-c()
準(zhǔn)備好循環(huán)體-可使用嵌套的for循環(huán)，完成完整的計算但不重復(fù)
注意第二個for循環(huán)的值為 i: nrow(data)，這一點(diǎn)很巧妙for (i in 1:nrow(data)){
for (r in i:nrow(data)){
g1=rownames(data)[i]
g2=rownames(data)[r]
c_r=cor(as.numeric(data[i,]),as.numeric(data[r,]),method="pearson")
p=cor.test(as.numeric(data[i,]),as.numeric(data[r,]),method ="pearson")[[3]]
##保存每一步的數(shù)據(jù)，而不可直接以空向量作為每一步運(yùn)行的結(jié)果
gene_name1=c(gene_name1,g1)
gene_name2=c(gene_name2,g2)
cor_r=c(cor_r,c_r)
pvalue=c(pvalue,p)
}
}
輸出為數(shù)據(jù)框，導(dǎo)出結(jié)果 data_cor<-data.frame(gene_name1,gene_name2,cor_r,pvalue)
head(data_cor)
## gene_name1 gene_name2 cor_r pvalue
## 1 gene1 gene1 1.0000000 7.395571e-32
## 2 gene1 gene2 0.4436884 3.781395e-01
## 3 gene1 gene3 0.2553650 6.252788e-01
## 4 gene1 gene4 0.2900609 5.771108e-01
## 5 gene1 gene5 -0.3356649 5.154125e-01
## 6 gene1 gene6 -0.6095215 1.989414e-01
dim(data_cor)
## [1] 210 4
|
