Unsupervised Learning and NGS¶

class(Golub_Merge)
head(Golub_Merge)
summary(Golub_Merge)

#Fortunately, bioconductor has a function to extract the expression data:

Golub.expression.data=exprs(Golub_Merge)
head(Golub.expression.data)
dim(Golub.expression.data)

# Now, we compute the distance matrix:

Golub.distance.matrix=dist(t(Golub.expression.data))

# We used 't' to transpose the matrix. This is just because dist works on rows and we want distance on columns.

as.matrix(Golub.distance.matrix)

dim(as.matrix(Golub.distance.matrix))

# cmdscale takes the distance matrix and produces a 72x2 matrix,
#    while preserving the distances between the 72 subjects.

MDS=cmdscale(Golub.distance.matrix)
as.matrix(dist(MDS))[,1:3]

# Now, we just plot the result.

plot(MDS,col=c(2,4)[ pData(Golub_Merge)[["ALL.AML"]]],pch=19,type="n",xlab="",ylab="")
text(MDS,as.character(pData(Golub_Merge)[["ALL.AML"]]),col=c(2,4)[pData(Golub_Merge)[["ALL.AML"]]])

# We observe the clustering of AML patients together and ALL patients together.

# Simulate noisy HM
# load the necessary libraries

library(pheatmap)
library(genefilter)

# n = number of members of each group, m = number of genes, alpha = significance level
# The data simulation is almost exactly the same as the noise simulation from the supervised learning lab

simulate.noise.heatmap=function(n,m,alpha)
  {
    # Simulate Expression Matrix

    EXPRS=matrix(rnorm(2*n*m),m,2*n) # Just randomly chosen, normally distributed values (N(0,1))

    rownames(EXPRS)=paste("Gene",1:m,sep="")
    colnames(EXPRS)=paste("patient id",1:(2*n),sep="")

    # Assign first n patients to group 1, second n patients to group 2
    grp=factor(rep(0:1,c(n,n)))

    # Get the two sample t-statistics
    # run rowttests on the expression values and extract the p-values for each gene
    pvals=rowttests(EXPRS, grp)$p.value


    #Choose the genes that have p-values smaller than significance level
    topgenes=which(pvals<alpha)
    EXPRS=EXPRS[topgenes,]

    # Just annotate which group each patient is in
    annodat=data.frame(Condition=ifelse(grp==0,"N","Y"),row.names=colnames(EXPRS))


    # draw heatmap
    pheatmap(EXPRS,
             border_color =NA,
             show_rownames = FALSE,
             show_colnames=FALSE,
             annotation_col=annodat,
             color=colorRampPalette(c("red3", "black", "green3"))(50),
             annotation_colors=list(Condition=c(Y="blue",N="yellow")))

    return(length(topgenes))  # return the number of topgenes found
}

set.seed(765)

# run with sample of 20000 genes and significance 0.005
aa=simulate.noise.heatmap(20,20000,0.005)
print(aa)

#Now simulate 40000 genes with significance 0.0025
aa=simulate.noise.heatmap(20,40000,0.0025)
print(aa)

# And now a sample of only 6 patients, with 20000 genes and significance level of 0.005
aa=simulate.noise.heatmap(3,20000,0.005)
print(aa)