HWr2

See also the current and the the previous lecture.

• Do either tasks A,B,C (beginners) or B,C,D (more advanced). You can also do all four for bonus credit.
• In your protocol write used R commands with brief comments on your approach.
• Submit required plots with filenames as specified.
• For each task also include results as required and a short discussion commenting the results/plots you have obtained. Is the value of interest increasing or decreasing with some parameter? Are the results as expected or surprising?
• Outline of protocol is in /tasks/r2/protocol.txt

Task A: sign test

• Consider a situation in which players played n games, out of which a fraction of q were won by A (the example in the lecture corresponds to q=0.6 and n=10)
• Compute a table of p-values for n=10,20,...,90,100 and for q=0.6, 0.7, 0.8, 0.9
• Plot the table using matplot (n is x-axis, one line for each value of q)
• Submit the plot in sign.png
• Discuss the values you have seen in the plot / table

Outline of the code:

# create vector rows with values 10,20,...,100
rows=(1:10)*10
# create vector columns with required values of q
columns=c(0.6, 0.7, 0.8, 0.9)
# create empty matrix of pvalues 
pvalues = matrix(0,length(rows),length(columns))
# TODO: fill in matrix pvalues using binom.test

# set names of rows and columns
rownames(pvalues)=rows
colnames(pvalues)=columns
# careful: pvalues[10,] is now 10th row, i.e. value for n=100, 
#          pvalues["10",] is the first row, i.e. value for n=10

# check that for n=10 and q=0.6 you get p-value 0.3769531
pvalues["10","0.6"]

# create x-axis matrix (as in HWr1, part D)
x=matrix(rep(rows,length(columns)),nrow=length(rows))
# matplot command
png("sign.png")
matplot(x,pvalues,type="l",col=c(1:length(columns)),lty=1)
legend("topright",legend=columns,col=c(1:length(columns)),lty=1)
dev.off()

Task B: Welch's t-test on microarray data

Read the microarray data, transform it to log scale, then work with table a:

input=read.table("/tasks/r2/acids.tsv", header=TRUE, row.names=1)
a = log(input)

Columns 1,2,3 are reference, columns 4,5,6 acetic acid, 7,8,9 benzoate, 10,11,12 propionate, and 13,14,15 sorbate

Write a function my.test which will take as arguments table a and 2 lists of columns (e.g. 1:3 and 4:6) and will run for each row of the table Welch's t-test of the first set of columns versus the second set. It will return the resulting vector of p-values, one for each gene.

• For example by calling pa <- my.test(a, 1:3, 4:6) we will compute p-values for differences between reference and acetic acid (computation may take some time)
• The first 5 values of pa should be

> pa[1:5]
[1] 0.94898907 0.07179619 0.24797684 0.48204100 0.23177496

• Run the test for all four acids
• Report how many genes were significant with p-value at most 0.01 for each acid
  • See Vector arithmetics in HWr1
  • You can count TRUE items in a vector of booleans by sum, e.g. sum(TRUE,FALSE,TRUE) is 2
• Report how many genes are significant for both acetic and benzoate acids simultaneously (logical and is written as &).

Task C: multiple testing correction

Run the following snippet of code, which works on the vector of p-values pa obtained for acetate in task B

# adjusts vectors of p-vales from tasks B for using Bonferroni correction
pa.adjusted = p.adjust(pa, method ="bonferroni")
# add this adjusted vector to frame a
a <-  cbind(a, pa.adjusted)
# create permutation ordered by pa.adjusted
oa = order(pa.adjusted)
# select from table five rows with the lowest pa.adjusted (using vector oa)
# and display columns containing reference, acetate and adjusted p-value
a[oa[1:5],c(1:6,16)]

You should get an output like this:

            ref1     ref2     ref3  acetate1   acetate2  acetate3 pa.adjusted
SUL1    7.581312 7.394985 7.412040 2.1633230 2.05412373 1.9169226 0.004793318
YMR244W 2.985682 2.975530 3.054001 0.3364722 0.33647224 0.1823216 0.188582576
DIP5    6.943991 7.147795 7.296955 0.6931472 0.09531018 0.5306283 0.253995075
YLR460C 5.620401 5.801212 5.502482 3.2425924 3.48431229 3.3843903 0.307639012
HXT4    2.821379 3.049273 2.772589 7.7893717 8.24446541 8.3041980 0.573813502

Do the same procedure for benzoate p-values and report the result (in your table, report both p-values and expression levels for bezoate, not acetate). Comment the results for both acids.

Task D: volcano plot, test on data generated from null hypothesis

Draw a volcano plot for the acetate data

• The x-axis of this plot is the difference between the mean of reference and the mean of acetate. You can compute row means of a matrix by rowMeans.
• The y-axis is -log10 of the p-value (use the p-values before multiple testing correction)
• You can quickly see the genes that have low p-values (high on y-axis) and also big difference in the mean expression between the two conditions (far from 0 on x-axis). You can also see if acetate increases or decreases the expression of these genes.

Now create a simulated dataset sharing some features of the real data but observing the null hypothesis that the mean of reference and acetate are the same for each gene

• Compute vector m of means for columns 1:6 from matrix a
• Compute vectors sr and sa of standard deviations for reference columns and for acetate columns respectively. You can compute standard deviation for each row of a matrix by apply(some.matrix, 1, sd)
• For each i in 1:6398, create three samples from the normal distribution with mean m[i] and standard deviation sr[i] and three samples with mean m[i] and deviation sa[i] (use the rnorm function)
• On the resulting matrix apply Welch's t-test and draw the volcano plot.
• How many random genes have p-value at most 0.01? Is it roughly what we would expect under the null hypothesis?

Draw a histogram of p-values from the real data (reference vs acetate) and from the random data (use function hist). Describe how they differ. Is it what you would expect?

Submit plots volcano-real.png, volcano-random.png, hist-real.png, hist-random.png (real for real expression data and random for generated data)