Hi everyone,
Here is the thing, we have recently submitted a paper about tumor prognosis biomarker, and a reviewer asked us to validate our findings in a given microarray dataset on http://www.ebi.ac.uk/arrayexpress.
After I downloaded the Raw data (.CEL files) and opened it with R (which I have never used before), I can only find manuals of full-scale expression analysis which took hours on my computer. After that, I have to pick out that specific gene's expression, integrate it with the survival data, and go for the validations, which is quite time-consuming.
I am wondering if there is a way to only extract the expression data of that specific gene(e.g. gene name 284065_at in HG U133 plus 2) from the Raw data files(in this case .CELs) and get the results I need much faster.
Meanwhile, does anyone knows how can I combine that expression data directly with survival data, (or better, direct data processing online without using R on my computer)? Or if you know a way to improve my work flow(look below for the codes I used), please share with us :)
P.S. the R code I am using is from this site: http://jura.wi.mit.edu/bio/education/bioinfo2007/arrays/array_exercises_1R.html
first installed Bioconductor, then:
1.Use the pull-down menu (File >> Change dir [Windows]) to go to the directory where you put the raw data (CEL files).
2.Load the "library" that contains the Affymetrix microarray code we'll want to use with the command
library(affy)
3.Read the CEL files (first command below) and then summarize and normalize with MAS5 (second command below). This could take a few minutes.
affy.data = ReadAffy()
eset.mas5 = mas5(affy.data)
4.The variable 'eset.mas5' contains normalized expression values for all probesets, along with other information. Let's continue by getting the expression matrix (probesets/genes in rows, chips in columns).
exprSet.nologs = exprs(eset.mas5)
# List the column (chip) names
colnames(exprSet.nologs)
5.exprSet = log(exprSet.nologs, 2) (not quite sure what this does)
6.To print out our expression matrix (as with most data), we can use a command like
write.table(exprSet, file="RESULTS_mas5_matrix.txt", quote=F, sep="\t")
to get a tab-delimited file that we could view in Excel or a text editor.
7.While we're doing Affymetrix-specific preprocessing, let's calculate an Absent/Present call for each probeset.
# Run the Affy A/P call algorithm on the CEL files we processed above
data.mas5calls = mas5calls(affy.data)
# Get the actual A/P calls
data.mas5calls.calls = exprs(data.mas5calls)
# Print the calls as a matrix
write.table(data.mas5calls.calls, file="Su_mas5calls.txt", quote=F, sep="\t")
PPS: these microarray datasets are really amazing, wish I understand more of the tools to utilize them~
Try this
# Here I am showing you an example of breast cancer data analysis that is on GSE27447
##go to R by just clicking R
R
#install the core bioconductor packages, if not already installed
source("http://bioconductor.org/biocLite.R")
biocLite()
# install additional bioconductor libraries, if not already installed
biocLite("GEOquery")
biocLite("affy")
biocLite("gcrma")
biocLite("hugene10stv1cdf")
biocLite("hugene10stv1probe")
biocLite("hugene10stprobeset.db")
biocLite("hugene10sttranscriptcluster.db")
#Load the necessary libraries
library(GEOquery)
library(affy)
library(gcrma)
library(hugene10stv1cdf)
library(hugene10stv1probe)
library(hugene10stprobeset.db)
library(hugene10sttranscriptcluster.db)
## see which directory you are in
getwd()
#Set working directory for download
setwd("/your_directory_path")
#Download the CEL file package for this dataset (by GSE - Geo series id)
getGEOSuppFiles("GSE27447")
#Unpack the CEL files
setwd("./your_directory_path/GSE27447")
untar("GSE27447_RAW.tar", exdir="data")
cels = list.files("data/", pattern = "CEL")
sapply(paste("data", cels, sep="/"), gunzip)
cels = list.files("data/", pattern = "CEL")
setwd("./your_directory_path/GSE27447/data")
raw.data=ReadAffy(verbose=TRUE, filenames=cels, cdfname="hugene10stv1") #From bioconductor
#perform RMA normalization (I would normally use GCRMA but it did not work with this chip)
data.rma.norm=rma(raw.data)
#Get the important stuff out of the data - the expression estimates for each array
rma=exprs(data.rma.norm)
#Format values to 5 decimal places
rma=format(rma, digits=5)
#Map probe sets to gene symbols or other annotations
#To see all available mappings for this platform
ls("package:hugene10stprobeset.db") #Annotations at the exon probeset level
ls("package:hugene10sttranscriptcluster.db") #Annotations at the transcript-cluster level (more gene-centric view)
#Extract probe ids, entrez symbols, and entrez ids
probes=row.names(rma)
Symbols = unlist(mget(probes, hugene10sttranscriptclusterSYMBOL, ifnotfound=NA))
Entrez_IDs = unlist(mget(probes, hugene10sttranscriptclusterENTREZID, ifnotfound=NA))
#Combine gene annotations with raw data
rma=cbind(probes,Symbols,Entrez_IDs,rma)
#Write RMA-normalized, mapped data to file
write.table(rma, file = "rma.txt", quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)
# here you would find the normalized expression of your gene
# if possible make boxplot to see if data is properly normalized or not
boxplot(rma, las=2)
## tell me how it looks
Hi Manvendra, apologies that I didn't reply earlier. I have tried your method, but seems it also gets the expression of the whole array (every gene tested) first, the same as the method I used.
I am looking for a method to skip the whole array-data extraction step (took very long time on my PC), and extract only the target gene I needed. Still, I am very grateful for your answer!
- Badges
- Science topic
- Similar topics
- Biological Science
- Genetics
- Gene Expression