Standardization

x <- matrix(1:21, ncol=7)

By row (gene)
std.x.by.row <- t(scale(t(x), scale=T))

By column (array)

std.x.by.col <- scale(x, scale=T)

Check the batch effect by date

date2col <- function(date.list)
{
clr.template = c("red", "orange", "yellow", "green", "cyan", "blue", "purple")
num.dates <- length(date.list)
clr.list <- vector()
clr.list[1] <- "red"
c.index <- 0
for (i in 2:num.dates) {
if(date.list[i] == date.list[i-1]) {
clr.list[i] = clr.list[i-1]
} else {
c.index <- c.index+1
clr.list[i] = clr.template[c.index %% 7+1]
}
}
return(clr.list)
}



setwd("E:\\CEL\\GastricCancer\\AU\\PM_data_new\\Gastric_Affy_files\\Tumors")
data <- read.table(file="AU_GC70.rma.txt", header=T, row.names=1)
data.ctrl <- data[54614:54675, ]
library("gplots")


#data <- t(scale(t(data.ctrl), scale=T)) #standardized by row(gene)
#data[data < -3] <- -3
#data[data > 3] <- 3

data <- sweep(data.ctrl, 1, apply(data.ctrl, 1, median)) #just median centered

my.color <- c("8/4/2004","8/4/2004","11/18/2004","11/18/2004","11/25/2004","11/25/2004","11/25/2004","11/25/2004","11/25/2004","11/25/2004","11/25/2004","11/26/2004","11/26/2004","11/26/2004","11/26/2004","11/26/2004","12/2/2004","12/2/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","12/3/2004","1/14/2005","1/14/2005","1/14/2005","1/14/2005","1/14/2005","1/14/2005","2/17/2005","2/17/2005","2/17/2005","2/17/2005","2/17/2005","2/17/2005","2/17/2005","2/25/2005","2/25/2005","3/4/2005","3/4/2005","3/4/2005","3/18/2005","3/18/2005","3/23/2005","4/8/2005","4/8/2005","4/8/2005","4/8/2005","4/8/2005","4/8/2005","4/28/2005","4/28/2005","4/28/2005","4/28/2005","4/29/2005","4/29/2005","4/29/2005","4/29/2005","5/19/2005","5/24/2005","5/24/2005","6/22/2005","6/22/2005","6/22/2005","6/22/2005","6/22/2005")
#my.color <- rep("black",dim(data)[2])
my.color <- date2col(my.color)
hm<-heatmap.2(as.matrix(data), col=greenred(75), scale="none", dendrogram="none", Rowv= T, Colv=F, ColSideColors=my.color, key=TRUE, symkey=FALSE, density.info="none",trace="none", cexRow=0.75,cexCol=0.75)
pdf(file = "Batch_in_CtrlGenes.pdf", width=10, height=10)
#pdf(file = "Batch_in_CtrlGenes.pdf")
hm<-heatmap.2(as.matrix(data), col=greenred(75), scale="none", dendrogram="none", Rowv= T, Colv=F, ColSideColors=my.color, key=TRUE, symkey=FALSE, density.info="none",trace="none", cexRow=0.75,cexCol=0.75)
dev.off()

GSAA

High, mid, low activity, e.g. p53
mid vs low -->(GSEA)  ES
high vs low --> High ES??

Supposed we have 200 cell lines, obtain the expression before drug treatment, and obtain GI50 for M drugs.

(NCI60, too small?)

	Drug1	Drug2	Drug3	... ...	DrugM
GeneSet1	Corr(1,1) = GSAA1vsGI50 for drug 1 across all cell lines	Corr(1,2) = GSAA1vsGI50 for drug2 across all cell lines
GeneSet2
GeneSet3
GeneSet4
... ...
GeneSetN

Median Centered

data.log10.MeanCentered <- sweep(data.log10, 2, apply(data.log10, 2, median))

Enrichment test in cancer subtypes

	Invasive	Proliferative	Metabolic
Hypermethylated	4110	1536	988
Hypomethylated	625	1263	284

Q1: Is the Invasive subtype enriched with Hypermethylated CpGs?

q <- 4110 # number of success (white balls) drawn

m <- 4110+1536+988 # number of success (white) in the urn

n <- 625+1263+284 # number of fail (black) in the urn

k <- 4110+625 # number of balls drawn out

p.value <- phyper(q-1,m,n,k, lower.tail=F)

# why q-1, because we want to have Pr(X>=x) instead of Pr(X>x)

# ref: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724271/ , search "phyper"

p.value

= 2.975069e-162

YES

Conclusion: The Invasive subtype is significantly enriched with Hypermethylated CpGs.

Q2: Is the Invasive subtype enriched with Hypomethylated CpGs?

q <- 625 # number of success (white balls) drawn

m <- 625+1263+284 # number of fail (black) in the urn

n <- 4110+1536+988 # number of success (white) in the urn

k <- 4110+625 # number of balls drawn out

p.value <- phyper(q-1,m,n,k, lower.tail=F)

p.value

=1

No

Q3: Is the Proliferative subtype enriched with Hypermethylated CpGs?

q <- 1536

m <- 4110+1536+988

n <- 625+1263+284

k <- 1263+1536

p.value <- phyper(q-1,m,n,k, lower.tail=F)

p.value

 = 1

No

Q4: Is the Proliferative subtype enriched with Hypomethylated CpGs?

q <- 1263

m <- 625+1263+284

n <- 4110+1536+988

k <- 1263+1536

p.value <- phyper(q-1,m,n,k, lower.tail=F)

p.value

= 3.765696e-193

YES

Conclusion: The Proliferative subtype is significantly enriched with Hypomethylated CpGs.

Q5: Is the Metabolic subtype enriched with Hypermethylated CpGs?

q <- 988 # number of success (white balls) drawn

m <- 4110+1536+988 # number of success (white) in the urn

n <- 625+1263+284 # number of fail (black) in the urn

k <- 988+284 # number of balls drawn out

p.value <- phyper(q-1,m,n,k, lower.tail=F)

p.value

= 0.01919936

YES

Conclusion: The Metabolic subtype is significantly enriched with Hypermethylated CpGs.

Q6: Is the Metabolic subtype enriched with Hypomethylated CpGs?

q <- 284 #

m <- 625+1263+284

n <- 4110+1536+988

k <- 988+284

p.value <- phyper(q-1,m,n,k, lower.tail=F)

p.value

= 0.9839128

No

See also:
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2724271/
An integrative genomics approach identifies Hypoxia Inducible Factor-1 (HIF-1)-target genes that form the core response to hypoxia

Median survival time

smed <- function(x) {

        ox <- capture.output(print(x))
        n <- length(ox)
        tmp <- t(sapply(ox[4:n],        
                function(l) strsplit(l, split=' +')[[1]]))
        nres <- strsplit(ox[3],split=' +')[[1]][2:6]
        res <- matrix(as.numeric(tmp[,2:6]), ncol=5,
                dimnames=list(tmp[,1], nres))
        res

}

# example:
#library(survival)
#fit1 <- survfit(Surv(time, status) ~ 1, data=aml)
#sf1 <- smed(fit1)
#sf1
#fit <- survfit(Surv(time, status) ~ x, data=aml)
#sf <- smed(fit)
#sf

http://tolstoy.newcastle.edu.au/R/help/06/05/26713.html

By Xiaochun Li

Change the $PATH

[leiz@136-25 ~]$ echo $PATH
/usr/lib64/qt-3.3/bin:/usr/kerberos/bin:/usr/local/bin:/bin:/usr/bin:/usr/NX/bin:/home/leiz/bin:/home/leiz/NGS/bowtie:/home/leiz/NGS/bowtie/indexes:/home/leiz/NGS/tophat:/home/leiz/NGS/samtools:/home/leiz/NGS/cufflinks

[leiz@136-25 ~]$ more .bash_profile
# .bash_profile

# Get the aliases and functions
if [ -f ~/.bashrc ]; then
        . ~/.bashrc
fi

# User specific environment and startup programs

PATH=$PATH:$HOME/bin:/home/leiz/NGS/bowtie:/home/leiz/NGS/bowtie/indexes:/home/leiz/NGS/tophat:/home/leiz/NGS/samtools:$HOME/NGS/cufflinks


export PATH

vi
source .bash_profile