Глубина покрытия в 1кб окнах

Question

Я хотел бы изобразить среднюю глубину охвата по моему геному, с хромосомами, выстроенными в порядке возрастания. Я вычислил глубину покрытия на позицию для моего генома, используя samtools. Я хотел бы создать график (который использует окна 1 КБ), как показано на рисунке 7: http://www.g3journal.org/content/ggg/6/8/2421/F7.large.jpg?width=800&height=600&carousel=1

Пример кадра данных:

Chr   locus depth
chr1    1   20  
chr1    2   24  
chr1    3   26  
chr2    1   53  
chr2    2   71  
chr2    3   74  
chr3    1   29  
chr3    2   36  
chr3    3   39  

Нужно ли менять формат кадра данных, чтобы разрешить непрерывную нумерацию для переменной V2? Есть ли способ усреднить каждые 1000 строк, и построить 1kb окон? И как бы я занялся черчением?

ОБНОВЛЕНИЕ РЕДАКТИРОВАНИЯ: Я смог создать новый набор данных как скользящее среднее из неперекрывающихся окон размером 1 Кбайт, используя эту запись: Покрытие генома как скользящее окно , и я сделал V2 непрерывным, т. Е. (1: 9 вместо 1,2,3 , 1,2,3,1,2,3)

library(reshape) # to rename columns
library(data.table) # to make sliding window dataframe
library(zoo) # to apply rolling function for sliding window

#genome coverage as sliding window
Xdepth.average<-setDT(Xdepth)[, .(
  window.start = rollapply(locus, width=1000, by=1000, FUN=min, align="left", partial=TRUE),
  window.end = rollapply(locus, width=1000, by=1000, FUN=max, align="left", partial=TRUE),
  coverage = rollapply(coverage, width=1000, by=1000, FUN=mean, align="left", partial=TRUE)
), .(Chr)]

И на сюжет

library(ggplot2)
Xdepth.average.plot <- ggplot(Xdepth.average, aes(x=window.end, y=coverage, colour=Chr)) + 
  geom_point(shape = 20, size = 1) +
  scale_x_continuous(name="Genomic Position (bp)", limits=c(0, 12071326), labels = scales::scientific) +
  scale_y_continuous(name="Average Coverage Depth", limits=c(0, 200))

Мне не повезло с использованием facet_grid, поэтому я добавил опорные линии, используя geom_vline(xintercept = c(). Смотрите ответ, который я разместил ниже для получения дополнительной информации / кодов, а также ссылок на графики. Теперь мне просто нужно поработать над маркировкой ...

Answer 1

Чтобы ответить на часть вопроса, вы пытались добавить + facet_grid(~ Chr) к вашему графику?(или + facet_grid(~ V2) в зависимости от имен ваших переменных)

Я не вижу вашего сообщения об ошибке, если я использую ваши примеры данных.Сообщение часто появляется, когда вы пытаетесь принять log(0), поэтому вы можете добавить псевдосчет log(x + 1), принять преобразование sqrt или asinh (последнее, если вы используете отрицательные значения).Что касается примеров данных, это хороший этикет, чтобы опубликовать пример данных в формате, который могут быть скопированы другими пользователями для проверки вашей проблемы, например:

depth <- data.frame(
  Chr = paste0("chr", c(1, 1, 1, 2, 2, 2, 3, 3, 3)),
  locus = c(1, 2, 3, 1, 2, 3, 1, 2, 3),
  depth = c(20, 24, 26, 53, 71, 74, 29, 36, 39)
)

Чтобы обратиться к части биоинформатики, выВозможно, вы захотите взглянуть на пакет GenomicRanges bioconductor: он имеет функцию tileGenome() для создания бинов, и вы можете использовать findOverlaps() со своими данными и бинами.Если у вас есть эти перекрытия, вы можете split() свои данные, основываясь на том, что они перекрывают, и рассчитать среднее покрытие для каждого разделения.

Имейте в виду, что вам может потребоваться некоторое время, чтобы ознакомиться сGRanges структура объекта и получение данных в этом (или GPos) формате.GRanges объекты похожи на файлы постели с геномными интервалами, в то время как GPos объекты похожи на точные, одиночные нуклеотидные координаты.

Однако вы уверены, что не хотите, чтобы число считываний на одну ячейку вместо среднего покрытия?Следует иметь в виду, что покрытие немного смещено против длинных операций чтения.

В качестве решения, не относящегося к R, вы также можете использовать bamCoverage в наборе deeptools с размером, скажем, 1000bp.

РЕДАКТИРОВАТЬ: воспроизводимый пример для построения

library(ggplot2, verbose = F, quietly = T)
suppressPackageStartupMessages(library(GenomicRanges))

# Setting up some dummy data
seqinfo <- rtracklayer::SeqinfoForUCSCGenome("hg19")
seqinfo <- keepStandardChromosomes(seqinfo)
granges <- tileGenome(seqinfo, tilewidth = 1e6, cut.last.tile.in.chrom = T)
granges$y <- rnorm(length(granges))

# Convert to dataframe
df <- as.data.frame(granges)

# The plotting
ggplot(df, aes(x = (start + end)/2, y = y)) +
  geom_point() +
  facet_grid(~ seqnames, scales = "free_x", space = "free_x") +
  scale_x_continuous(expand = c(0,0)) +
  theme(aspect.ratio = NULL,
        panel.spacing = unit(0, "mm"))

^{Создано в 2019-04-22 с помощью пакета Представить (v0.2.1)}

Answer 2

Играя больше с программой, я смог создать новый набор данных в виде скользящего среднего значения не перекрывающихся окон размером 1 Кбайт, используя этот пост: Покрытие генома в виде скользящего окна , которое не заняло много времени или не затягивало много памяти.

library(reshape) # to rename columns
library(data.table) # to make sliding window dataframe
library(zoo) # to apply rolling function for sliding window
library(ggplot2)

 #upload data to dataframe, rename headers, make locus continuous, create subsets
depth <- read.table("sorted.depth", sep="\t", header=F)
depth<-rename(depth,c(V1="Chr", V2="locus", V3="coverageX", V3="coverageY")
depth$locus <- 1:12157105
Xdepth<-subset(depth, select = c("Chr", "locus","coverageX"))

#genome coverage as sliding window
Xdepth.average<-setDT(Xdepth)[, .(
  window.start = rollapply(locus, width=1000, by=1000, FUN=min, align="left", partial=TRUE),
  window.end = rollapply(locus, width=1000, by=1000, FUN=max, align="left", partial=TRUE),
  coverage = rollapply(coverage, width=1000, by=1000, FUN=mean, align="left", partial=TRUE)
), .(Chr)]

Для построения нового набора данных:

#plot sliding window by end position and coverage
Xdepth.average.plot <- ggplot(Xdepth.average, aes(x=window.end, y=coverage, colour=Chr)) + 
  geom_point(shape = 20, size = 1) +
  scale_x_continuous(name="Genomic Position (bp)", limits=c(0, 12071326), labels = scales::scientific) +
  scale_y_continuous(name="Average Coverage Depth", limits=c(0, 250))

Затем я попытался добавить facet_grid(. ~ Chr) для разделения по хромосоме, но каждая панель отстоит далеко друг от друга и повторяет всю ось вместо того, чтобы быть непрерывной.

Обновление: я пробовал различные настройки с scales = "free_x" и space = "free_x". Наиболее близким было удаление пределов из scale_x_continuous() и использование scales = "free_x" и space = "free_x" с facet_grid, но ширина панели по-прежнему не пропорциональна размеру хромосомы, а ось x очень шаткая. Для сравнения я вручную добавил контрольные линии, используя geom_vline(xintercept = c() между хромосомами (ожидаемый результат).

Идеальное разделение и ось X без надписей на панели с использованием

Xdepth.average.plot +
  geom_vline(xintercept = c(230218, 1043402, 1360022, 2891955, 3468829, 3738990, 4829930, 5392573, 5832461, 6578212, 7245028, 8323205, 9247636, 10031969, 11123260, 12071326, 12157105))

График с контрольными линиями

Снятие лимита с scale_x_continuous() и использование facet_grid

Xdepth.average.plot5 <- ggplot(Xdepth.average, aes(x=window.end, y=coverage, colour=Chr)) + 
  geom_point(shape = 20, size = 1) +
  scale_x_continuous(name="Genomic Position (bp)", labels = scales::scientific, breaks = 
                       c(0, 2000000, 4000000, 6000000, 8000000, 10000000, 12000000)) +
  scale_y_continuous(name="Average Coverage Depth", limits=c(0, 200), breaks = c(0, 50, 100, 150, 200, 300, 400, 500)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(legend.position="none")
X.p5 <- Xdepth.average.plot5 + facet_grid(. ~ Chr, labeller=chr_labeller, space="free_x", scales = "free_x")+
  theme(panel.spacing.x = grid::unit(0, "cm"))
X.p5

График с гранями и без ограничений по оси X