I'm working on a BioStats141 homework investigating the correlation between gene expression and cancer. There are a set of 22283 genes we have data for from 107 individuals. So each of the 107 individuals expresses each of the 22283 to varying degrees. If an individual has a higher expression of Gene# than others, and they smoke, are they more likely to get cancer?
I've been trying to piece together this question since 9am, and now I'm starting to be overwhelmed by all the data and different correlations I'm finding. I'm not quite sure which data is the important ones that I need to correlate, so I'm going to take a break, have lunch, print out the question on paper and see if I can't annotate and get some new insight.
I am a master of looping now, haha! Show you the code later.
# Homework 7
# Stats 141 / Bio 141
# Now we extract the information we need for our analysis:
expressions = exprs(my.eset) # extracts the matrix of gene expressions
# one row for each "feature" (genes)
# one column for each sample (individual)
# 107 individuals, 22283 rows "feature" (genes)
phenotypes = pData(my.eset) # extracts a data.frame containing the
# descriptions of the samples (individuals)
# 107 individuals (rows), 7 descriptors (colums)
annotation = fData(my.eset) # extracts the annotation information for the platform
# used (Affymetrix Human Genome U133A Array)
# 22283 features (genes), 21 components
# phenotypes now contains a description of each sample. We will perform
# an ANOVA analysis, gene by gene (that is, one test for each row of
# the "expressions" matrix), based on the following categorical
# variables (extracted from the data frame "phenotypes"):
# "tissue" "individual" "disease.state" "gender"
> levels(phenotypes$tissue)
[1] "normal" "tumor"
> levels(phenotypes$individual)
[1] "current smoker" "former smoker" "never smoker"
> levels(phenotypes$disease.state)
[1] "stage I" "stage II" "stage III" "stage IV"
> levels(phenotypes$gender)
[1] "female" "male"
# We are interested in genes that show an interaction bewteen tissue
# and individual.
# So, we want to perform a test in which the null hypothesis is that
# the interaction term is 0, versus the alternative hypothesis that
# it is not, in the presence of the other predictors.
# Question: Why is this interesting? What distinguishes genes that
# are chosen by having a small p-value in this test?
A: In this test, genes with small p-values are affected by the interaction
term of smoking status and tissue type. These genes are likely to be
cancer-causing genes and should be further investigated.
###### ANOVA Test ######
# Separate the categorical columns from Phenotypes we're studying
tissue.type = data.frame("Tissue" = phenotypes$tissue)
smoking.status = data.frame("Status" = phenotypes$individual)
disease.state = data.frame("State" = phenotypes$disease.state)
gender = data.frame("Gender" = phenotypes$gender)
my.phenotypes <- cbind(tissue.type,smoking.status,disease.state,gender)
# Since the axis between Phenotypes and Expressions are switched, we have
# to rebind genes so that each gene is now a column instead of a row.
# We continue to consider each row representing an individual, as in
# Phenotypes.
for (i in 1:nrow(expressions)){
gene <- expressions[i,]
# for each gene...
Gene = data.frame("Gene" = gene)
# put into data frame
data <- cbind(my.phenotypes,Gene)
# bind to phenotype columns so that we can perform lm()
fit1 <- lm(gene ~ Tissue + Status,data=data)
# no interaction term
fit2 <- lm(gene ~ Tissue + Status + Tissue:Status,data=data)
# with interaction term
my.anova <- anova(fit1,fit2)
pp[i] <- my.anova$Pr[2]
# significance values of interaction term
}
###### Benjamini-Hochberg Procedure ######
# Choose an FDR level so that the Benjamini-Hochberg procedure yields
# approximately 100 discoveries. Report this FDR.
# Put pp into a data frame to perform Hochberg
# Attach annotations to view later (when things are jumbled)
Hochberg = data.frame("PValue" = pp)
Annot = data.frame("Annotations" = annotation$Gene.title)
Hochberg2 <- cbind(Hochberg,Annot)
# Each gene now has a pvalue and its annotation only
# Sort our p-values
sorted.Hochberg <- Hochberg2[order(Hochberg2$PValue),]
# Attach i-values
i = c(1:22283)
final.Hochberg <- cbind(sorted.Hochberg,i)
# Save annotations order now that pvalues are sorted
final.Annot = data.frame("Annotations" = final.Hochberg$Annotations)
# State our hypotheses for the Benjamini-Hochberg procedure
# Ho: fit1 = fit2; interaction term is zero
# Ha: fit1 does not equal fit2; interaction term not zero
# keep trying alpha values until we get 100 discoveries
# Create a loop for the Benjamini-Hochberg procedure
alpha = .1911
m = 22283
for (i in 1:m){
TF[i] <- final.Hochberg[i,1] < (alpha*final.Hochberg[i,3])/m
}
sum(as.numeric(Hochberg.test$value))
[1] 100
# FDR alpha value of 0.1911 or higher will give us 100 discoveries
# TF contains a set of TRUE FALSE values, where TRUE indicates that
# the Ho test was rejected. FALSE values are those that fail to
# reject Ho.
# Put TF into a data frame which we can bind to sorted.data
# Bind annotations before sorting TRUE/FALSE values
Hochberg.test = data.frame("value" = TF)
Hochberg.test2 <- cbind(Hochberg.test,final.Annot)
# Sort the Hochberg test to separate TRUE/FALSE and their annotations
sorted.Hochbergtest <- Hochberg.test2[order(Hochberg.test2$value),]
# Since FALSE values = 0, the last 100 are TRUE values and are the selected
# genes
sorted.Hochbergtest[22183:22283,]
# Briefly examine the annotation of the selected genes by using
# annotation$Gene.title
# and comment on the results.
A: These are the genes that seem to be affected by the interaction term,
under the Benjamini-Hochberg test with FDR alpha = 0.1911. Their names
do not immediately mean anything to me, but we could use this test to
investigate the correlation of these genes to smoking and tumors.
---------------------------------------------------------
Ok, so the analysis is terrible. I had no idea what I was talking about. But the codes are cool. And long.
No comments:
Post a Comment