Showing posts with label Stat530-hw2. Show all posts
Showing posts with label Stat530-hw2. Show all posts

2007-09-26

The main dish is coming? Q1.1 ~ 4... XD

Q1.1) fold change:
easy thought, easy Done! Sorting data by AVE(schz)/AVE(ctrl) and then pick up the first 200 genes. But how?


>aa<-array(round(rnorm(24),2),dim=c(6,4))  # randomly generate number based on normal dist
>aa
[,1] [,2] [,3] [,4]
[1,] -0.53 -1.26 0.87 1.92
[2,] 1.00 0.26 -0.40 -0.67
[3,] 1.39 -0.16 -0.11 -0.06
[4,] -0.51 0.35 0.28 0.07
[5,] -1.59 1.80 0.52 -0.63
[6,] 0.03 -0.60 0.61 -0.17
>order(aa[ ,2]) # sort by column 2
>order(aa[3, ]) # sort by row 5
>order(aa[3, ], decreasing=TRUE) # this one is the default
>order(aa[3, ], decreasing=FALSE)
>aa[order(aa[3,]),]
>aa[order(aa[3,])[1:2],3:5] #aa[order(aa[3,])[row],col]
As you can see the order() returning the order of row or column as a vector, and that's why it's a beautiful thing for us when selecting genes by kinda order.

Heatmap is so pretty and vivid.. XD

Don't just try heatmap..There are more options such as heatmap.2 from gplots package, and heatmap_2 from Heatplus package.
Both heatmap.2 and heatmap_2 can show the bar of color code for arbitrary unit (color key ); in addition to that, heatmap_2 can repress the dendrogram for rows (but unable to suppress the clustering) and show different colors for dendrogram of columns based on their clusters.
Don't worry about the data format! They "eat" the same format of data, meaning what heatmap can load is good for heatmap.2 & heatmap_2 too.
Some miscellaneous setup just for pretty look
1) color pattern
Color palette can be easily assigned by using col=ColorPattern.
Here are some internal/pre-defined color palettes:
heat.colors(n), topo.colors(n), redgreen(n) and cm.colors(n)
Change (n) to set up how many levels you want in the color gradient
rev(): reverse the color pattern, for example rev(redgreen(100)) will get a gradient from green to red instead of the one from red to green.

2) color-stripe for specific group of sample
scol=c(...) is a vector with the length equaling to the number of your column.
Just assign the color you want in th vector, like c("red","red","skyblue","skyblue") for a 4-coulmned dataset.
3) scaling
scaling="row", "col", "none"; although keeping this option ON usually makes the heatmap pretty, I thought turning it OFF can help us have an ideal how many genes is in relatively low expression levels.

How to feed the program ?

The first challenge met is what kind of function I should use to feed the program.
Checked the data file, even with .xls; actually it's a tab separated text file, and thus I decided to use read.table() function.

Read.table has very flexible parameters allowed to set up for reading in titles/column names, change of separators(^t, ^p, "." and so on). However, R also provides several subtypes of read.table, which with mild difference of default.

One of them, read.delim(), fits to my needs. It reads in the first row as the column names (check it by using colnames()), recognizes "tab" as separator (type ?read.delim to see detail).

So far so good, right? Don't be happy too early..
Since our purpose is to visualize the data by heatmap, I'm so eager to see it first.. (what?! just put the fold change or something else on the back strove). I'm so naive that it's not a big deal to throw the data into the heatmap(). Here comes the problem... Orz

>t<-read.delim("sch.dat") >t
>colnames(t)
>dim(t)
>heatmap(t[2:7072,6:29],col=rev(heat.colors(100)))
Error in heatmap(t[2:7072, 6:29], col = rev(heat.colors(100))) :
'x' must be a numeric matrix
#since the column(1:5) are not numeric, they are probe/gene information
#It's so weird, the range I select t[2:7072, 6:29] should be ALL numeric.
#NOW learn a lesson what you saw is NOT what you think!!!
> t[1,7:8]
C2.A.. C3.A..
1 18291.58 19402.7
> is.numeric(t[1,7:8]) # check if it is numeric.. XD XD XD it's NOT.. Orz.. no wonder..
[1] FALSE
>NewSet <- array(rnorm(7071*25),dim=c(7071,25)) # put some number to array (just for sure it's numeric!!!)
>for(i in 1:7071){for(j in 1:24){NewSet[i,j]=t[i,5+j]}} # transfer data from original array to new one
>NewSet[1,2:3] # look almost identical to t[1,7:8]
[1] 18291.58 19402.70
>is.numeric(NewSet[1,2:3])
[1] TRUE
>rownames(NewSet)
NULL
>colnames(NewSet)
NULL
#now the data matrix has no row/col names, which is not good for us to know which is which!
#assign the names to them
>rownames(NewSet)=t[,1] # row=genes
>colnames(NewSet)<-c("c1","c2","c3","c4","c5","c6","c7","c8","c9","c10","c11","c12","p1","p2","p3","p4","p5","p6","p7","p8","p9","p10","p11","p12","") ## The reason I would like to assign the row/col names is the heatmap function can automatically pick up the names and show on the plot
## OK! now the data is ready for heatmap !!!!?