# В файле counts5.txt лежат результаты RNA-seq: количество сиквенсов (ридов)
# детектированных для каждого гена (строчки) в каждом образце (столбец)
# всего 8 образцов: 4 контроля и 4 после обработки дигидрогена монооксидом (ДГМО)

setwd('/Users/tepliakova/Desktop/univer/r_course')
data = read.table(file='counts5.txt')
expr.ratio = data / lapply(data, sum)

# задача: найти гены значимо меняющие экспрессию (долю ридов приходящихся
# на них в образце) после обработки ДГМО, проконтролировав долю ложных предсказаний (FDR)

before = expr.ratio[c('c1', 'c2', 'c3', 'c4')]
after = expr.ratio[c('t1', 't2', 't3', 't4')]

p.values = numeric(dim(before)[1])
for (i in 1:dim(before)[1]) {
  p.values[i] = wilcox.test(as.numeric(before[i,]), as.numeric(after[i,]))$p.value
}

genes = 1:(dim(before)[1])
as.numeric(table(p.values < 0.05)['TRUE'])
fdr = 0.05 * dim(before)[1] / as.numeric(table(p.values < 0.05)['TRUE'])
# number of genes = 1337
# fdr = 0.748

# Используйте т-тест. Так как не известно насколько нормально распределение
# сделайте 100 пермутаций (перемешайте столбцы в таблице), и посчитайте сколько генов
# будет значимо в пермутации на том же уровне значимости, что вы выбрали в предыдущем пункте

permute = function(data,n=100){
  num.genes.valuable = numeric(n)
  for(i in 1:n){
    data = sample(data)
    p.values.t = numeric(dim(data)[1])
    for (j in 1:dim(data)[1]) {
      p.values.t[j] = t.test(as.numeric(data[j, 1:4]), as.numeric(data[j, -(1:4)]))$p.value
    }
    num.genes.valuable[i] = as.numeric(table(p.values.t < 0.05)['TRUE'])
  }
  num.genes.valuable
}

num.genes = permute(expr.ratio)

# Формат сдачи: число значимых генов, FDR. Весь необходимый код.