Презентация: http://makarich.fbb.msu.ru/artemov/R/R_6.pdf

Дифференциальная экспрессия генов

Загрузка данных http://makarich.fbb.msu.ru/artemov/R/counts.tab

counts=read.table("http://makarich.fbb.msu.ru/artemov/R/counts.tab")
head(counts)

Для поиска дифференциально экспрессирующихся генов потребуется пакет DESeq из набора пакетов bioconductoR.

http://bioconductor.org/packages/release/bioc/html/DESeq.html

Пакеты bioconductoR устанавливаются не так, как пакеты из стандартного репозитария:

source("http://bioconductor.org/biocLite.R")
#устанавливаем пакет
biocLite("DESeq")
#подключаем библиотеку
library(DESeq)

Создадим объект, в котором хранятся данные для теста дифференциальной экспрессии

condition = factor( c( "untreated", "untreated", "treated", "treated" ) )
cds = newCountDataSet( counts, condition )

Нормализация

Оценим поправочные коэффициенты для учета разного количества отсеквенированных в каждом образце ридов. Обратите внимание на синтаксис: данная операция что-то добавляет к объекту и возвращает новый измененный объект (в данном случае - с вычисленными коэффициентами). Таким образом, чтобы сохранить новый объект с тем же именем, используется запись вида x=f(x)

cds = estimateSizeFactors( cds )
#коэффициенты
sizeFactors( cds )
#как посмотреть на таблицу значений экспресси после нормировки?
countsNorm= counts( cds, normalized=TRUE )
head(countsNorm)

Тест

#оценим дисперсию для каждого гена: дисперсия оценивается как функция от средней экспрессии.
cds= estimateDispersions( cds)
#собственно, тест. Сравниваем группы untreated и treated. На выходе - таблица
res = nbinomTest( cds, "untreated", "treated" )
res[res$id=="FBgn0026562",]
#можно упорядочить по убыванию
xr=res[order(res$padj),]
#уберем строки в таблице, где есть NA (неопределенные значения)
xr=na.omit(xr)
length(which(xr$padj<0.05))

Полезный способ визуализации:

plotMA(res)

Поправка на множественное тестирование

Что такой p adjusted или FDR? Если среди 10000 генов нашлись 500 генов с p-value<0.05 - это не удивительное событие. Даже если бы никакого эффекта не было, то доля 0.05 всех генов имела бы p-value<0.05, таким образом, мы бы увидели 10000*0.05=500 генов, которые ошибочно посчитали бы значимыми.

Выход: ужесточать (уменьшать) порог на p-value или умножать p-value на какие-то коэффициенты (увеличивать). Второй подход - и есть поправка на множественное тестирование.

http://makarich.fbb.msu.ru/artemov/R/MultipleTesting.pdf

res$pval - (не скорректированные) значения p-value. Попробуем их вручную скорректировать

#(очень жесткая) поправка Бонферрони
padj_Bonferroni=p.adjust(res$pval, method="bonferroni")
#поправка Benjamini-Hochberg, или FDR. Обысно используется она
padj_FDR=p.adjust(res$pval, method="fdr")
#её ли использует DESeq?
all( is.na(res$pval) | res$pval==padj_FDR )

Работа с геномными разметками

Главный инструмент - библиотека GenomicRanges, позволяет работать с геномными разметками и эффективно искать пересечения отрезков.

source("http://bioconductor.org/biocLite.R")
biocLite("GenomicRanges")
library("GenomicRanges")

Загрузим разметку экзонов дрозофилы

biocLite("biomaRt")
library(biomaRt) 
#http://www.ensembl.org/biomart/martview
#какую базу данных использовать
ensembl <- useMart("metazoa_mart_22",dataset="dmelanogaster_eg_gene")
#какие поля вытащить
fields <- c("chromosome_name","strand","ensembl_gene_id","ensembl_exon_id","start_position","end_position","exon_chrom_start","exon_chrom_end")
genes <- getBM( fields, mart=ensembl, filters="chromosome_name", values=c("4"))
head( genes )

Создадим GRanges из полученной таблицы

annot <- GRanges(
seqnames = Rle("chr4"),#хромосома
ranges=IRanges(
start=genes$exon_chrom_start, 
end=genes$exon_chrom_end
), #объект IRanges
strand = Rle(genes$strand), 
exon=genes$ensembl_exon_id, 
gene=genes$ensembl_gene_id
#все остальные поля
)

Чтение bam-файла и вычисление покрытия экзонов ридами

biocLite("Rsamtools")
library("Rsamtools")
aln_all <- readGappedAlignments("http://makarich.fbb.msu.ru/artemov/R/untreated1_chr4.bam")
head(aln_all)

overlaps=countOverlaps( annot, aln_all )
head(overlaps)

Дальше сложим получившиеся значения для всех экзонов каждого гена. В данном случае tapply выбирает группы значений overlaps, которым соответствую одинаковые значения annot$gene и применяет к группам функцию sum (то есть суммирует покрытие всех экзонов гена).

geneexpr=tapply(overlaps, list( annot$gene ), FUN=sum )
head(geneexpr)

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

Аннотация SNP

http://makarich.fbb.msu.ru/artemov/R/SNP_chr22.vcf

Загрузим файл с найденными полиморфизмами (формат VCF). Для этого потребуется пакет VariantAnnotation

source("http://bioconductor.org/biocLite.R")
biocLite("VariantAnnotation")
library("VariantAnnotation")

vcf = readVcf("/home/artemov/R/SNP_chr22.vcf", "hg19")

Загрузим bed-файл c экзонами человека.

Для этого понадобится пакет rtracklayer

source("http://bioconductor.org/biocLite.R")
biocLite("rtracklayer")
library("rtracklayer")

exons=import.bed("http://makarich.fbb.msu.ru/artemov/R/exons_chr22.bed")

При помощи функции countOverlaps посчитаем для каждого снипа, сколько раз он пересекается с экзонами. Нас будет интересовать, 0 или более 0 раз снип пересекся с экзоном.

overlapsExon=countOverlaps(vcf, exons)
sum(overlapsExon>0)
length(overlapsExon)

Выберем нужные полиморфизмы и сохраним их в файл

vcf_take=vcf[overlapsExon>0,]
dim(vcf_take)
writeVcf(vcf_take, "SNP_in_exons.vcf")

Kodomo

Пользователь