Анализ микрочиповых данных в R/ Bioconductor

Анализ микрочиповых данных в R/Bioconductor Александр Ишкин 20.12.2010

R и Bioconductor. Установка. • Установка R: • www.r-project.org • Установка Bioconductor • http://www.bioconductor.org/ • На данный момент полный набор – более 400 пакетов «на все случаи жизни». • Код для установки базового набора пакетов из R: • Этого хватит для большинства задач по анализу экспрессионных микрочипов. • source("http://bioconductor.org/biocLite.R") • biocLite()

Микрочипы. Различные типы. • Практически все типы чипов можно анализировать при помощи Bioconductor: • Экспрессия генов • Affymetrix GeneChip • Illumina BeadChip • Прочие (в том числе кДНК-микрочипы) • Вариабельность генома • Копийность ДНК (aCGH) • SNP-чипы • микроРНК-чипы • Метилирование • Экзонные чипы • Выстилающие (tiling)чипы • Проблема в том, чтобы найти подходящие пакеты • http://www.bioconductor.org/help/bioc-views/

Экспрессионные чипы. Последовательность действий

Нормализация • Огромное количество способов • Самый прижившийся – квантильная нормализация

Трансформация данных • Суть в том, чтобы убрать зависимость ст. отклонения от величины сигнала – она нарушает предположения многих тестов • В основном используется логарифмирование

Структуры данных • Для всех пакетов Bioconductor, имеющих дело с чипами, основная структура данных – класс ExpressionSet • x <- exprs(dataset) – получить матрицу с сигналами экспрессии • y <- pData(phenoData(dataset)) – получить аннотацию образцов • z <- featureNames(dataset) – получить названия зондов ExpressionSet + всякие метаданные sampleNames exprs Сигналы phenoData Аннотация образцов featureData Аннотация зондов featureNames

Тестовый массив данных • Взят из базы данных NCBI GEO; ID: GSE2737 • Чип: Affymetrix U95A. 22283 зонда (“probe sets”) • Экспрессия генов в псориатических бляшках • 4 образца – пораженная кожа, 4 – непораженная кожа, 3 – кожа здоровых людей (эти 3 на другом чипе)

Загрузка данных Affy • Загружаем «низкоуровневые» данные – CEL файлы. • CDF – определяют состав “probe set” • В исходных CDF 20% последовательностей зондов не выравниваются с транскриптами • Custom CDF - http://brainarray.mbni.med.umich.edu/Brainarray/Database/CustomCDF/CDF_download.asp ##1. Load Affymetrixdata library(affy) celfiles <- list.files("affy_raw", full.names=T) affy.raw <- ReadAffy(filenames = celfiles) ##Set alternative CDF affy.raw@cdfName <- "HGU95A_Hs_ENTREZG"

QC - Affy • Специфические контрольные параметры для Affy library(simpleaffy) affy.raw<- ReadAffy(filenames = celfiles) affy.qc<- qc(affy.raw) plot(affy.qc)

QC - Affy library(arrayQualityMetrics) affy.raw<- ReadAffy(filenames = celfiles) arrayQualityMetrics(affy.raw, outdir = “qc”) ##this wrote pdf report to “qc” directory • Более общие критерии – насколько образцы похожи друг на друга • Корреляция образцов • Сходство распределений сигналов

Препроцессинг - Affy • RMA делает следующие вещи: • Вычитание фонового шума • Квантильная нормализация сигналов • Логарифмирование (трансформация) • Суммирование сигналов (несколько сигналов “probes”  1 сигнал “probe set”, в нашем случае – сигнал с гена.) • В общем случае нормализацию и трансформацию надо делать отдельно • Для Affymetrix есть много таких алгоритмов (еще один – GCRMA). Кто из них самый-самый – вопрос открытый. ##Set alternative CDF affy.raw@cdfName <- "HGU95A_Hs_ENTREZG” ##Run summarization algorithm – RMA affy.data <- rma(affy.raw)

Препроцессинг- Affy • После RMA имеем массив данных, с которым можно работать • Запись в файл write.exprs(affy.data, "data/affy_rma.txt")

Анализ дифференциальной экспрессии • Для начала необходимо определить группы образцов, которые мы будем сравнивать. • В нашем случае можно задать их прямо в коде. Обычно лучше скомпоновать отдельный файл с аннотацией. • Если задача состоит в поиске дифференциальной экспрессии, то это самый ответственный момент. ##4. Compose phenotypic data disease.state <- rep(c("Inactive", "Active", 4) pdata <- data.frame(DiseaseState=disease.state, row.names=sampleNames(affy.data))

Анализ дифференциальной экспрессии

Анализ дифференциальной экспрессии • Код для вышеописанной последовательности действий ##5. Perform differential expression analysis ##Create design library(limma) design <- model.matrix(~ -1+pdata$DiseaseState) colnames(design) <- c("A", "I") ##Create contrasts contrast <- makeContrasts(A-I, levels=design) ##Fit model fit1 <- lmFit(affy.data, design) fit2 <- contrasts.fit(fit1,contrast) ##Apply moderated t-test with empirical Bayes correction fit3 <- eBayes(fit2) ##Get differentially expressed genes de.all <- topTable(fit3, coef=1, adjust="BH", number = nrow(affy.data)) de.genes <- de.all[de.all$adj.P.Val < 0.05, c("ID", "adj.P.Val")]

Анализ дифференциальной экспрессии • Результаты: • Запись в файл ##Write ##Convert IDs to Entrez Gene IDs de.genes$ID <- sub("_at", "", de.genes$ID) write.table(de.genes, "processed/affy_DE genes.txt", quote=F, row.names=T)

Кластеризация • Задача: исследовать структуру массива данных, выявить группы схожих генов/образцов • Иерархическая кластеризация ##6. Clustering. affy.dist <- dist(t(exprs(affy.data)), method="euclidean") affy.clust <- hclust(affy.dist, method="average") plot(affy.clust)

Кластеризация • Задача: исследовать структуру массива данных, выявить группы схожих генов/образцов • Principal Component Analysis (PCA) ##PCA affy.pca <- prcomp(t(exprs(affy.data))) col=rep(c("blue", "red"),4) label <- pdata$DiseaseState plot(affy.pca$x[,1:2], pch=19, main ="GSE2737 PCA", col=col) text(affy.pca$x[,1], affy.pca$x[,2]-1, labels=label, col=col)

Загрузка данных Illumina • Packages: beadarray, lumi • Достаточно простая последовательность действий, если данные представлены в «родном» форматеBeadStudio ##Import Illumina data library(lumi) illumina.file <- list.files("illumina_raw", full.names=T) ##Read raw data illumina.raw <- lumiR(illumina.file) ##Perform QC illumina.qc <- LumiQ(illumina.raw) ##Perform variance stabilization (analogous to log2) illumina.vst <- lumiT(illumina.raw) ##Perform quantile normalization illumina.data <- lumiN(illumina.vst) write.exprs(illumina.data, "data/illumina.txt")

Обработка данных по copy number variation • Суть: для каждого зонда посчитать общий относительный сигнал и выделить участки с последовательными координированными изменениями экспрессии • ПрепроцессингaCGHи SNP-чипов: crlmm, oligo, DNAcopy, aroma.affymetrix, snpChip • Сегментация: DNAcopy, BioHMM, RankCopy

Дополнительные возможности BioC • Работа с данными next-gen sequencing • Biostrings, ShortRead, BSGenome • chipseq(анализ данных ChIP-Seq) • Genominator, Deseq(анализ данных RNA-Seq) • Загрузка данных и метаданных из баз NCBI через R • GEOMetadb, GEOquery(GEO – микрочипы, RNA-seq) • SRAdb(NCBI Sequence Reads Archive - секвенирование) • Аннотация • GOstats, biomaRt, microarray annotation packages (hgu95a) • Интерфейс UCSC Genome Browser • rtracklayer • Анализ биологических сетей • BioNet, RCytoscape

Дополнительные возможности R • Взаимодействие с другими языками программирования • rJava; • rcpp • 4 • Взаимодействие с базами данных • ROracle, RSQLite…; • Параллелизация • R/parallel; multicore • Работа с большими массивами данных • http://www.slideshare.net/bytemining/r-hpc • Bigmemory, biglm

Обучающие материалы по R • http://en.wikibooks.org/wiki/R_Programming • Introduction to R by Longhow Lam - http://www.splusbook.com/RIntro/RCourse.pdf • http://www.statmethods.net/index.html - Quick-R, сборник материалов по применению статистических методов в R

Редакторы для R • IDE • Emacs + ESS (www.ess.r-project.org/) • Eclipse plugin (www.walware.de/goto/statet) • Текстовые редакторы • Tinn-R (www.sciviews.org/Tinn-R/) • Notepad++ (www.notepad-plus-plus.org/) • Во всех этих программах • Подсветка синтаксиса R • Запуск кода в консоль R

Анализ микрочиповых данных в R/ Bioconductor