R與bioconductor--ExpressionSet SummarizedExperiment GEOquery biomaRt S4-Classes S4-Methods

博主自學了coursera上來自約翰霍普金斯大學<使用Bioconductor分析基因組科學資料>，很不錯，推薦給大家 ExpressionSet Overview ExpressionSet -Expression Matrix -Phenotype data Feature data eSet：僅僅是有多個Expression Matrix


ExpressionSet
library(ALL) data(ALL) experimentData(ALL) exprs(ALL)[1:4,1:4] head(sampleNames(ALL)) head(featureNames(ALL)) head(pData(ALL)) ALL$sex ALL[1:10,1:5] featureData(ALL)#包含關於基因的資訊，但是經常沒有 ids <- featureNames(ALL)[1:5] library(hgu95av2.db)#之前資料集裡面寫了晶片平台hgu95av2.db as.list(hgu95av2ENTREZID[ids]) phenoData(ALL)#不如使用pData(ALL),phenoData(ALL)內容更多些 names(pData(ALL))#有時相當於varLabels(ALL),有時varLabels(ALL)更詳細

SummarizedExperiment
#SummarizedExperiment library(airway) data("airway") airway colData(airway)#pData(ALL) airway$cell colnames(airway) head(rownames(airway)) assayNames(airway)#要想獲得表達矩陣,使用assay accessor,用assayNames()獲得全部表達矩陣名 assay(airway,"counts")[1:4,1:4] length(rowRanges(airway)) rowRanges(airway)#SummarizedExperiment特別之處在於每行每列都有關聯的GRanges start(airway) gr = GRanges("1",ranges = IRanges(start = 1,end = 10^7)) subsetByOverlaps(airway,gr) #'正是有著關聯GRanges,可以取出染色體"某一地區內"基因的表達值

GEOquery
library(GEOquery) eList <- getGEO("GSE11675") length(eList) eData = eList[[1]] eData names(pData(eData)) eList2 = getGEOSuppFiles("GSE11675")#下載原始tar包 eList2


biomaRt
library(biomaRt) head(listMarts()) mart <- useMart("ensembl") mart head(listDatasets(mart)) ensemble <- useDataset("hsapiens_gene_ensembl",mart) values <- c("202763_at","209310_s_at","207500_at") getBM(attributes = c("ensembl_gene_id","affy_hg_u133_plus_2"),       filters = "affy_hg_u133_plus_2",values = values,mart = ensemble) attributes <- listAttributes(ensemble)#可以查詢到的條目 nrow(attributes)#可以把一個物種的基因轉換到另一個物種的同源基因 head(attributes) filters <- listFilters(ensemble)#可以查詢到的條目 nrow(filters)#可以把一個物種的基因轉換到另一個物種的同源基因 head(filters) attributePages(ensemble)#attributes儲存在一個一個page中,可以用這個減小搜尋範圍 attributes <- listAttributes(ensemble,page = "feature_page") nrow(attributes)
R S4 Classes
library(ALL) library(GenomicRanges) #'S3對象就是像一個list,list中每個對象都有各自的name #'而S4對象定義了每個class應該是有些什麼東西 data("ALL") ALL class(ALL) isS4(ALL) class?ExpressionSet#查看一個class的簡介 ?"ExpressionSet-class"#查看一個class的簡介 #'list的規則:首字母大寫 #'構造方法 ExpressionSet() getClass("ExpressionSet")#Slots插槽,就是這個class由哪兒些小class構成 ALL@annotation annotation(ALL) #class升級了,定義改變了,用updateObject OLD_OBJECT = updateObject(OLD_OBJECT) validObject(ALL)#檢測對象是否正確,是否符合class的定義

R S4 Methods
library(GenomicRanges) GenomicRanges::as.data.frame#S4方法 base::as.data.frame#S3方法 showMethods("as.data.frame")#可以看見,X類型不同,後續選用的程式碼也不同 #查看傳入某一特定類型，對應的相關程式碼 getMethod("as.data.frame","GenomicRanges") getMethod("as.data.frame",signature(x="GenomicRanges")) #查看傳入某一特定類型，對應的協助文檔 method?"as.data.frame,DataFrame" method?"as.data.frame,GenomicRanges" ?"as.data.frame,DataFrame-method" ?"as.data.frame,GenomicRanges-method" showMethods("findOverlaps") getMethod("findOverlaps",signature(query = "Ranges",subject = "Ranges")) ?"findOverlaps,Ranges,Ranges-method" #'S4缺點:難以找到help文檔,難以直接看原始碼,難以debug #'但是最好S4寫一個package,方便管理
最後是完整程式碼片段

#ExpressionSetlibrary(ALL)data(ALL)experimentData(ALL)exprs(ALL)[1:4,1:4]head(sampleNames(ALL))head(featureNames(ALL))head(pData(ALL))ALL$sexALL[1:10,1:5]featureData(ALL)#包含關於基因的資訊，但是經常沒有ids <- featureNames(ALL)[1:5]library(hgu95av2.db)#之前資料集裡面寫了晶片平台hgu95av2.dbas.list(hgu95av2ENTREZID[ids])phenoData(ALL)#不如使用pData(ALL),phenoData(ALL)內容更多些names(pData(ALL))#有時相當於varLabels(ALL),有時varLabels(ALL)更詳細#SummarizedExperimentlibrary(airway)data("airway")airwaycolData(airway)#pData(ALL)airway$cellcolnames(airway)head(rownames(airway))assayNames(airway)#要想獲得表達矩陣,使用assay accessor,用assayNames()獲得全部表達矩陣名assay(airway,"counts")[1:4,1:4]length(rowRanges(airway))rowRanges(airway)#SummarizedExperiment特別之處在於每行每列都有關聯的GRangesstart(airway)gr = GRanges("1",ranges = IRanges(start = 1,end = 10^7))subsetByOverlaps(airway,gr)#'正是有著關聯GRanges,可以取出染色體"某一地區內"基因的表達值library(GEOquery)eList <- getGEO("GSE11675")length(eList)eData = eList[[1]]eDatanames(pData(eData))eList2 = getGEOSuppFiles("GSE11675")#下載原始tar包eList2library(biomaRt)head(listMarts())mart <- useMart("ensembl")marthead(listDatasets(mart))ensemble <- useDataset("hsapiens_gene_ensembl",mart)values <- c("202763_at","209310_s_at","207500_at")getBM(attributes = c("ensembl_gene_id","affy_hg_u133_plus_2"),      filters = "affy_hg_u133_plus_2",values = values,mart = ensemble)attributes <- listAttributes(ensemble)#可以查詢到的條目nrow(attributes)#可以把一個物種的基因轉換到另一個物種的同源基因head(attributes)filters <- listFilters(ensemble)#可以查詢到的條目nrow(filters)#可以把一個物種的基因轉換到另一個物種的同源基因head(filters)attributePages(ensemble)#attributes儲存在一個一個page中,可以用這個減小搜尋範圍attributes <- listAttributes(ensemble,page = "feature_page")nrow(attributes)library(ALL)library(GenomicRanges)#'S3對象就是像一個list,list中每個對象都有各自的name#'而S4對象定義了每個class應該是有些什麼東西data("ALL")ALLclass(ALL)isS4(ALL)class?ExpressionSet#查看一個class的簡介?"ExpressionSet-class"#查看一個class的簡介#'list的規則:首字母大寫#'構造方法ExpressionSet()getClass("ExpressionSet")#Slots插槽,就是這個class由哪兒些小class構成ALL@annotationannotation(ALL)#class升級了,定義改變了,用updateObjectOLD_OBJECT = updateObject(OLD_OBJECT)validObject(ALL)#檢測對象是否正確,是否符合class的定義library(GenomicRanges)GenomicRanges::as.data.frame#S4方法base::as.data.frame#S3方法showMethods("as.data.frame")#可以看見,X類型不同,後續選用的程式碼也不同#查看傳入某一特定類型，對應的相關程式碼getMethod("as.data.frame","GenomicRanges")getMethod("as.data.frame",signature(x="GenomicRanges"))#查看傳入某一特定類型，對應的協助文檔method?"as.data.frame,DataFrame"method?"as.data.frame,GenomicRanges"?"as.data.frame,DataFrame-method"?"as.data.frame,GenomicRanges-method"showMethods("findOverlaps")getMethod("findOverlaps",signature(query = "Ranges",subject = "Ranges"))?"findOverlaps,Ranges,Ranges-method"#'S4缺點:難以找到help文檔,難以直接看原始碼,難以debug#'但是最好S4寫一個package,方便管理