Лекции по биоинформатике: Анализ экспрессии. Демонстрация.

Слайды лекции: (pdf)

Тема

Данные экспрессии генов являются одними из важнейших объектов исследования современной биоинформатики. Цель данного курса — рассказать о данных экспрессии и основных методах анализа таких данных.

Некоторые аспекты темы я постараюсь проиллюстрировать практическими примерами на R. При желании и должном умении быстро печатать и следить за лектором и клавиатурой одновременно, приведенные примеры можно будет повторить "вживую". Прошу заметить, что это не обязательно — примеры с лекции будут выложены здесь для желающих повторить их самостоятельно, — но если вы все равно планируете сидеть с открытым лаптопом то лучше уж открывайте его по делу :).

Для выполнения демонстраций необходимо поставить себе систему R и несколько пакетов к ней. Лучше сделать это заранее, а то в лекционном зале на всех может не хватить интернета.

Подготовка

Скачайте R:
- Windows: (тут)
- Линукс: что-нибудь в духе apt-get install r-base.
- Другие платформы: (тут)
Запустите R. Убедитесь что консоль работает.
```
2+2
1:10
for (i in 1:10) { print(i) }
help()
?plot
example(plot)
```
В будущем (скорее всего не на этой лекции) вам может пригодиться вот эта бумажка.
Поставьте закачиваться несколько пакетов из набора "Биокондуктор" (они довольно объемистые и это займет время). Обратите внимание, что в идеале нижеприведенные строки можно просто скопировать и вставить (Ctrl+V) прямо в консоль.

source("http://www.bioconductor.org/biocLite.R")
biocLite("ygs98probe")
biocLite("ygs98.db")
biocLite("BSgenome.Scerevisiae.UCSC.sacCer1") # ~12 мегабайт
biocLite("Biostrings")
biocLite("ArrayExpress")
biocLite("affy")
biocLite("annotate")
biocLite("GO.db")                             # ~ 17 мегабайт
biocLite("GOstats")
biocLite("seqLogo")
install.packages("impute")  # Это уже не биокондуктор, но тоже пригодится
install.packages("ggplot2") # Это на лекции не понадобится, это просто реклама

NB: Если вы подсоединены к интернету через прокси сервер, проставьте системную переменную окружения: http_proxy=http://vash.proxy.host.ru:3128 и перезапустите R, иначе пакеты не будут скачиваться.
Пока все устанавливается, гляньте на список пакетов в Биокондукторе

Часть I: Мыс Технологий

Affymetrix Probes

library(ygs98probe)
ls("package:ygs98probe")
ygs98probe
?ygs98probe
as.data.frame(ygs98probe[1:20,])

library(ygs98.db)
ls("package:ygs98.db")
ygs98ORF[["10000_at"]]
ygs98CHR[["10000_at"]]
ygs98CHRLOC[["10000_at"]]
ygs98CHRLOCEND[["10000_at"]]

library(BSgenome.Scerevisiae.UCSC.sacCer1)
ls("package:BSgenome.Scerevisiae.UCSC.sacCer1")
Scerevisiae
Scerevisiae$chr12

probe = ygs98probe[1,1]
seq = Scerevisiae$chr12[790669:791046]

library(Biostrings)
matchPattern(probe, seq)
matchPattern(probe, reverseComplement(seq))

as.data.frame(ygs98probe[1,]) # 102

ArrayExpress

library(ArrayExpress)
library(affy)
affydata = ArrayExpress("E-GEOD-18037") # ~8 мегабайт
affydata

pm(affydata)[1:10,]
mm(affydata)[1:10,]

dim(pm(affydata))
dim(ygs98probe)

Качество данных микрочипа

# График зависимости PM от ММ
smp = sample(1:nrow(pm(affydata)), 10000)
plot(pm(affydata)[,1], mm(affydata)[,1], xlab="PM", ylab="MM")
plot(log(pm(affydata)[smp,1]), log(mm(affydata)[smp,1]), xlab="log(PM)", ylab="log(MM)")

# График зависимости интенсивности зонда от его GC-содержания
# [Взято из: openVignette("Biostrings") -> третий пункт]
all_probes = DNAStringSet(ygs98probe[,1])
alphabets = alphabetFrequency(all_probes, baseOnly=TRUE)
gc_content = ordered(alphabets[,"C"] + alphabets[,"G"])
intensity = pm(affydata)[,1]
boxplot(log(intensity) ~ gc_content, col=rainbow(nlevels(gc_content)))

# Зависимость сигнала от положения на чипе
probe_pos = floor(ygs98probe[,2,3]/100)
probe_pos$value = pm(affydata)[,1]
agg_probes = aggregate(probe_pos, by=list(probe_pos$x,probe_pos$y), FUN=mean)
mx = matrix(0, nrow=6, ncol=6)
for (i in 1:nrow(agg_probes)) mx[agg_probes$x[i]+1, agg_probes$y[i]+1] = agg_probes$value[i]
image(mx)

# Дополнительный материал (на лекции не затрагивается):
# biocLite("arrayQualityMetrics")  # Скачивает большую кучу всякого добра, осторожно.
# library(arrayQualityMetrics)
# openVignette("arrayQualityMetrics")  # -> см. статью номер 2

Финальная нормализация

# Можно так
library(affy)
expdata = mas5(affydata)

# А можно так
expdata = rma(affydata)

# Дальше читаем здесь:
# openVignette("affy")  # (1 - туториал, 2 - описание стандартных методов)

# Теперь у нас объект типа ExpressionSet
?ExpressionSet
expdata
exptable = get("exprs", assayData(expdata))

Часть II: а.о. Многомерного Анализа

Визуализации

# Качаем данные
expdata = read.table("http://genome-www.stanford.edu/cellcycle/data/rawdata/combined.txt",
        sep='\t', header = TRUE, row.names=1)
# NB: Те же данные можно получить в пакете yeastCC: biocLite("yeastCC")

expdata = expdata[,26:49]   # Оставляем только нужные нам микрочипы.
names(expdata)
dim(expdata)
expdata[1:5,1:5]

# Восполним отсутствующие значения
#  См: http://www-stat.stanford.edu/~hastie/Papers/missing.pdf
library(impute)
expdata = impute.knn(as.matrix(expdata))$data

# Возьмем выборку, чтобы не напрягать зазря процессор
set.seed(1)         # Чтобы у всех получился одинаковый результат
rsamp = sample(1:nrow(expdata), 1000)
expdata.sample = expdata[rsamp,]

# Порисуем..
image(t(expdata.sample))
heatmap(expdata.sample)
#
plot(expdata.sample[1,],type='l',xlab='Experiment',ylab='Measurement')
for (i in 2:50) lines(expdata.sample[i,],col=i)
#
plot(expdata.sample[,1], expdata.sample[,2], xlab=colnames(expdata)[1], ylab=colnames(expdata)[2])
#
pc = prcomp(expdata.sample, retx=TRUE)
plot(pc)
plot(pc$x[,1], pc$x[,2], xlab="PC1", ylab="PC2")

pc = prcomp(t(expdata.sample), retx=TRUE)
plot(pc$x[,1], pc$x[,2], xlab="PC1", ylab="PC2")
text(pc$x[,1],pc$x[,2],colnames(expdata),pos=3)

Кластеринг

c = hclust(dist(expdata.sample, "euclidean"), "complete")
plot(c)
c = kmeans(expdata.sample, 10)
plot(pc$x[,1], pc$x[,2], col=c$cluster, xlab="PC1", ylab="PC2", pch=20)

Финальный кластеринг

# NB: Здесь еще можно было бы сделать:
#   Выкинуть гены, экспрессия которых слабо меняется
#   Нормализовать экспрессию каждого гена (тогда k-means
#   по свойствам больше похож на кластеризацию по корреляции)
#   На лекции мы этого делать не будем.
# Наводка:
#   expdata.sd = apply(expdata,1,sd)
#   hist(expdata.sd)
#   expdata.filt = expdata[which(expdata.sd > 0.5),]
#
#   meanCenter = function(x) { (x - mean(x)) / sd(x) }
#   expdata.mc = t(apply(expdata.filt, 1, meanCenter))

set.seed(1)
c = kmeans(expdata, 12)
csize = table(c$cluster)
csize
par(mfrow=c(3,4))
for (i in 1:12) 
    plot(c$centers[i,], type='l', col=i, xlab="",ylab="",
         main=sprintf("Cluster %d (%d genes)",i,csize[i]))

# Возьмем для примера шестой кластер
my = which(c$cluster == 6)
expdata.my = expdata[my,]
# И посмотрим на него
par(mfrow=c(1,1))
plot(expdata.my[1,],type='l',xlab='Experiment',ylab='Measurement',axes=FALSE)
axis(1,at=1:24,labels=colnames(expdata.my))
axis(2)
for (i in 2:nrow(expdata.my)) lines(expdata.my[i,],col=i)

Часть III: земля Смысла

Генетическая онтология

gene_names = rownames(expdata.my)

library(org.Sc.sgd.db)
names(org.Sc.sgdGO[["YGR221C"]])

nrow(expdata)
go7117 = unique(org.Sc.sgdGO2ALLORFS[["GO:0007117"]])
length(go7117)
length(gene_names)
length(intersect(go7117, gene_names))

dhyper(1, 40, 6178-40, 97)
1-phyper(0, 40, 6178-40, 97)

GOstats

library(GOstats)

params = new("GOHyperGParams", 
   geneIds=gene_names,
   universeGeneIds=rownames(expdata), 
   annotation="org.Sc.sgd.db", 
   ontology="BP", 
   pvalueCutoff=0.01/1000, 
   testDirection="over", 
   conditional=TRUE)

hgTest = hyperGTest(params)
hgTest
head(summary(hgTest))

paste(gene_names, collapse=" ")  # Получившуюся строку можно вставить в 
                                 # g:Profiler (http://biit.cs.ut.ee/gprofiler)

Часть IV: республика Нуклеотидов

Промотеры

# Этот пакет содержит полный геном S.cerevisiae
library(BSgenome.Scerevisiae.UCSC.sacCer1)

# Функция, возвращающая промотер гена по его имени
# Промотер (здесь) = 500 нуклеотидов до начала гена (включая первый нуклеотид гена)
# Если данных нет, возвращает NULL
getPromoter = function(g) {
    chr = org.Sc.sgdCHR[[g]]                            # Узнаем на какой хромосоме ген
    if (is.null(chr)) return(NULL)
    g_start = org.Sc.sgdCHRLOC[[g]]                     # .. и в какой он позиции
    chrdata = Scerevisiae[[sprintf("chr%s", chr)]]      # Берем всю хромосому
    p_start = g_start - 499
    if (g_start < 0 & p_start < 0)                      # И вырезаем из нее нужный кусок
        return(complement(chrdata[abs(p_start):abs(g_start)]))
    else if (g_start > 0 & p_start > 0)
        return(chrdata[p_start:g_start])
    else              # Если оказалось что промотер залезает за начало или конец хромосомы
        return(NULL)  # Лень разбираться, скажем что не повезло
}

# Находим промотер для каждого элемента массива gene_names
ps = apply(as.matrix(gene_names),1,getPromoter)

# Сравнивать будем со множеством случайных генов (так проще)
set.seed(1)
random_genes = sample(rownames(expdata), 500)
ns = apply(as.matrix(random_genes),1,getPromoter)
names(ns) = random_genes # Потом пригодится

# Удалим пустые элементы
null_ps = unlist(lapply(ps,is.null))
ps[which(null_ps)] = NULL
null_ns = unlist(lapply(ns,is.null))
ns[which(null_ns)] = NULL

Поиск Мотивов

# Сгенерируем список всех 8-нуклеотидных строк, входящих
# в промотеры генов нашего кластера (мотивов-кандидатов)
strs = c()
for (i in 1:length(ps))
    for (j in (1:(length(ps[[i]])-7)))
        strs = c(strs, toString(ps[[i]][j:(j+7)]))
strs = unique(strs)

# Проиндексируем все мотивы-кандидаты, чтобы быстро их искать
pd = PDict(strs, tb.start=3, tb.width=3)

# Посчитаем в скольки последовательностях из кластера
# встречается каждый мотив-кандидат
countPos = rep(0, length(strs)) # Нули
for (i in 1:length(ps))
    countPos = countPos + 
        (countPDict(pd, ps[[i]], max.mismatch=2) > 0)

# Проделаем то же самое для набора случайных последовательностей
countNeg = rep(0, length(strs))
for (i in 1:length(ns))
    countNeg = countNeg + 
        (countPDict(pd, ns[[i]], max.mismatch=2) > 0)

# Теперь считаем для каждого мотива-кандидата P-value
pvalues = c()
for (i in 1:length(strs))
    pvalues[i] = 1-phyper(countPos[i]-1, length(ps), length(ns), countPos[i]+countNeg[i])

# Применяем поправку FDR
pv = p.adjust(pvalues, "fdr")

# Берем самые интересные мотивы
motifs = which(pv < 0.001)

Отрисовка лого

# Возьмем первый найденный мотив
m = strs[motifs[1]]

# Найдем все вхождения его с двумя ошибками
all_matches = as.character(matchPattern(m, ps[[1]], max.mismatch=2))
for (i in 2:length(ps))
    all_matches = c(all_matches, as.character(matchPattern(m, ps[[i]], max.mismatch=2)))

# Построим матрицу вероятностей по позициям (PWM)
cm = consensusMatrix(all_matches, as.prob=TRUE)[1:4,]

# ...и нарисуем картинку
library(seqLogo)
seqLogo(cm)

Заключение: Предсказывающие модели

Гипотетический пример

y = expdata[names(ns),1]
x1 = sapply(ns, function(s) { countPattern(m, s, max.mismatch=2) })
x2 = sapply(ns, function(s) { countPattern("GATTACAG", s, max.mismatch=2) })
boxplot(y~x1)
lm(y~x1+x2)

Чуть менее гипотетический пример
Не пробуйте его если у вашего компьютера менее ~3GB оперативной памяти. Вдобавок (если у вас всего, скажем, 4GB), имеет смысл перезапустить R чтобы освободить память от объектов сделанных в предыдущей сессии, и повыключать конкурирующие за память приложения.
Конечно же этот пример можно было бы реализовать намного более эффективным образом, но наша цель здесь — проиллюстрировать суть а не заниматься оптимизацией.

# Сгенерируем множество всех подстрок длиной 7:
pastex = function(x,y) { paste(x,y,sep="") }
all_motifs1 = c("A","C","T","G")
all_motifs2 = as.character(outer(all_motifs1, all_motifs1, pastex))
all_motifs3 = as.character(outer(all_motifs2, all_motifs1, pastex))
all_motifs4 = as.character(outer(all_motifs2, all_motifs2, pastex))
all_motifs7 = as.character(outer(all_motifs3, all_motifs4, pastex))

# Проиндексируем все последовательности
pd = PDict(all_motifs7, tb.start=3, tb.width=3)

# Возьмем в этот раз ВСЕ промотеры 
# (это займет некоторое время, не нервничайте сходите заварите чаю)
all_promoters = apply(as.matrix(rownames(expdata)),1,getPromoter)
names(all_promoters) = rownames(expdata)

# Удалим пустые элементы
null_promoters = unlist(lapply(all_promoters,is.null))
all_promoters[which(null_promoters)] = NULL

# Для каждого промотера проверим наличие каждого мотива (подождите еще немного)
getMotifVector = function(seq) { as.integer(countPDict(pd, seq, max.mismatch=1) > 0) }
motifVectors = lapply(all_promoters, getMotifVector)

# motifVectors - это список векторов.
# преобразуем его в матрицу
M = matrix(0, nrow=length(motifVectors), ncol=length(motifVectors[[1]]))
for (i in 1:nrow(M)) M[i,] = motifVectors[[i]]
rm(motifVectors)    # Освободим память

rownames(M) = names(all_promoters)
colnames(M) = all_motifs7

M[1:5,1:5]

g = expdata[rownames(M),1]
as.matrix(g[1:5], ncol=1)

# Посчитаем абсолютное значение корреляции мотивов с экспрессией
motifCorrs = rep(0,ncol(M))
for (i in 1:ncol(M)) motifCorrs[i] = abs(cor(M[,i], g))
names(motifCorrs) = colnames(M)
sort(motifCorrs, decreasing=TRUE)[1:5]

# Обратите внимание, что в разных экспериментах разный "главный мотив".
mts = rep("",ncol(expdata))
for (j in 1:ncol(expdata)) {
    motifCorrs = rep(0,ncol(M))
    for (i in 1:ncol(M)) motifCorrs[i] = abs(cor(M[,i], expdata[rownames(M),j]))
    names(motifCorrs) = colnames(M)
    mts[j] = names(sort(motifCorrs, decreasing=TRUE))[1]
}