Возьмем данные RNA-seq об экспрессии генов при раке молочной железы (breast cancer) и в соответствующих нормальных образцах (взятых из тех же пациентов). http://makarich.fbb.msu.ru/artemov/R/BRCA_paired_20samp.tab
(*) Для подготовки таблицы использовались данные TCGA (The Cancer Genome Atlas), https://tcga-data.nci.nih.gov/tcga/ Их удобнее загружать с ресурса, предоставляемого Broad institute: https://confluence.broadinstitute.org/display/GDAC/Download
Обычная загрузка таблицы (в каждой ячейке - counts - количество ридов на ген).
Единственная сложность - понять, какой образец опухоль, какой - нормальная ткань. Формат названий образцов для TCGA описан здесь: https://wiki.nci.nih.gov/display/TCGA/TCGA+Barcode
Четвертый раздел названия (напр., "01A") описывает тип образца. Для нас важно, что первая цифра 0 соответствует опухолям, 1 соответствует нормам. Эта первая цифра вытаскивается в строке, где функции strsplit (разделяет строку по определенному символу/подстроке) и substr (вытаскивает определенную подстроку, в данном случае - первый символ). Цифру мы потом переименовываем в N - normal или T - tumour. https://tcga-data.nci.nih.gov/datareports/codeTablesReport.htm?codeTable=sample%20type
mexp=read.table("BRCA_paired_20samp.tab")
cond=sapply(colnames(mexp), function(x){substr( strsplit(x, ".", fixed=T)[[1]][4],1,1) } )
cond[cond=="0"]="T"
cond[cond=="1"]="N"
condition = factor( cond )Установка пакета DESeq - не как обычно, т.к. пакет не из стандартного репозитория, а из bioconductoR:
source("http://bioconductor.org/biocLite.R")
biocLite("DESeq")Собственно, тест на дифференциальную экспрессию
library(DESeq) EX = newCountDataSet( mexp, condition ) EX = estimateSizeFactors( EX ) EX= estimateDispersions( EX ) DE = nbinomTest( EX, "N", "T" )
(*) Дополнительное задание: учесть парность образцов (пары опухоль - образец нормальной ткани взяты из одного пациента)
Упорядочите гены от наибольшего p-value до наименьшего, скопируйте список top дифференциально экспрессирующихся генов. Посмотрите, какие категории генов перепредставлены среди дифференциально экспрессируемых, например, с помощью сервиса DAVID.