setwd("f:\\geo\\geo晶元資料/")
load('gse35896_eset.rdata')
a=gset[[1]]
##取出第乙個元素賦值給乙個物件a
dat=exprs(a)
#a現在是乙個物件，取a這個物件通過看說明書知道要用exprs這個函式，該函式得到表達矩陣
#現在 得到的dat就是乙個表達矩陣，只不過基因的id是探針名
dim(dat)
#看一下dat這個矩陣的維度
dat[1:5,1:5]
##以下為gpl570的包
library(hgu133plus2.db)
ids=totable(hgu133plus2symbol)
head(ids)
colnames(ids)=c('probe_id','symbol')
length(unique(ids$symbol))
#[1] 18832個獨特的基因探針，意味著本來19825個裡面有一部分是重複的
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
#每個物件出現的個數
plot(table(sort(table(ids$symbol))))
#畫圖觀察
ids=ids[ids$symbol != '',]
ids=ids[ids$probe_id %in% rownames(dat),]
##%in%用於判斷是否匹配，然後取匹配的幾行，去掉無法匹配的資訊。
dat[1:5,1:5]
dat=dat[ids$probe_id,]
#取表達矩陣中可以與探針名匹配的那些，去掉無法匹配的表達資料，這時只剩下19825個探針及表達資訊，其餘已被剔除。
#ids新建median這一列，列名為median，同時對dat這個矩陣按行操作，取每一行的中位數，將結果給到median這一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = t),]
#對ids$symbol按照ids$median中位數從大到小排列的順序排序
##即先按symbol排序，相同的symbol再按照中位數從大到小排列，方便後續保留第乙個值。
##將對應的行賦值為乙個新的ids，這樣order()就相當於sort()
ids=ids[!duplicated(ids$symbol),]
#將symbol這一列取取出重複項，'!'為否，即取出不重複的項，去除重複的gene ，保留每個基因最大表達量結果.最後得到18832個基因。
dat=dat[ids$probe_id,]
#新的ids取出probe_id這一列，將dat按照取出的這一列中的每一行組成乙個新的dat
rownames(dat)=ids$symbol
#把ids的symbol這一列中的每一行給dat作為dat的行名
head(dat)
##至此我們就得到了該資料集的表達矩陣，最後將結果儲存。
dat[1:3,1:3]
which(rownames(dat)=="myc")
pd=pdata(a)
#通過檢視說明書知道取物件a裡的臨床資訊用pdata
view(pd)
## 檢視一下，挑選一些感興趣的臨床表型，這裡我們欲得到其分組title資訊。
library(stringr)
#執行乙個字元分割包
group_list=str_split(pd$source_name_ch1,' ',simplify = t)[,4]
#抽取title一列，按照空格分割，取第四個元素即control和vemurafenib
table(group_list)
#看一下兩個分組各有幾個
group_list<-as.data.frame(group_list)
group_list$id<-pd$geo_accession
rownames(group_list)<-group_list$id
which(group_list$group_list=="primary")
rownames(group_list[1:566,])
tumordat<-dat[,1:566]
colnames(tumordat)
#tumordat<-dat
###得到腫瘤表達矩陣,整理gsea檔案
all_df <- cbind(name = rownames(tumordat), description = "na", tumordat)
#all_df <- all_df[,-1]
dim(all_df)
all_df[1:3,1:3]
group <- c(rep("low", 31), rep("high", 31))
group <- paste(group, collapse = " ")
group <- c(paste(c(62, 2, 1), collapse = " "), "# low high", group)
write.table(file = "geo35896.txt", all_df, sep = "\t", col.names = t, row.names = f, quote = f)
write.table(file = "group35896.cls", group, col.names = f, row.names = f, quote = f)
###相關性分析
myc<-dat["myc",]
gapdh<-dat["gapdh",]
data<-rbind(myc,gapdh)
data<-t(data)
data<-as.data.frame(data)
class(data)
###library(ggplot2)
#求相關係數
r<-cor(data$myc,data$gapdh,method = "pearson")
#p值p<-cor.test(data$myc,data$gapdh,alternative="two.side")
##畫圖
tiff(filename = "35896_cor.tiff")
ggplot(data =data,aes(x=myc, y=gapdh)) +
geom_point(color="#d7191c")+
geom_smooth(method="lm",color="#1a9641") +
geom_text(aes(x=4.4, y=10,label=paste("r","=",signif(r,3),
seq="")),
color="black") +
geom_text(aes(x=4.4, y=9.7,label="p = 0.105"),color="black")
dev.off()
R語言相關性分析及步驟
記錄一下r語言學習過程，對於r的基礎就是基本資料型別 向量，矩陣，資料框，字串等等 庫的呼叫以及函式自定義，還需要多加學習！進入主題，今天主題是相關性分析 以下為 y c 170,175,180 定義向量 y1 c 20,25,30 y cor.test y,y1,method spearman 呼...
R語言相關性分析
相關性分析就是通過定量指標描述變數之間的強弱 直接或間接的聯絡。常見相關性指標 pearson相關係數是用於表示相關性大小的最常用指標，數值介於 1 1之間，越接近0相關性越低，越接近 1或1相關性越高。正負號表明相關方向，正號為正相關 負號為負相關。又稱為秩相關係數，利用兩變數的秩次大小來進行分析...
基於R進行相關性分析
一 相關性矩陣計算 1 載入資料 data read.csv 231 6057 2016 04 05 zx wd 2.csv header false 說明 csv格式的資料，header false 表示沒有標題，即資料從第一行開始。2 檢視匯入資料的前幾行，3 刪除資料的7，8列，都是0 4 計...