August 8th, 2013

R-Bioconductor - rtracklayer

Дисклемер - краткое описание функционала инструмента импорта/экспорта файлов геномных аннотаций.

rtracklayer используется как интерфейс между R, геномным браузером и геномными аннотациями. Главный функционал - ввод/вывод табулированных (в т.ч. сжатых) файлов геномных аннотаций: GTF/GFF, BED, WIG, BedGraph, BED15 и BigWig. Загружает файлы аннотаций в GRanges, способен различать разные версии GFF если указывать точный формат, WIG также хранит в виде GRanges.
Специфическая ошибка из-за чувствительности к разделителям в последних столбцах GFF показала баг в выводе UCSC Table Browser - лишние пробелы в концах строк.
Ввод/вывод сессий и треков браузера не интересен, гораздо важнее скриптовый веб-интерфейс к UCSC Table Browser: задаются базовые параметры сессии, и по имени треки и таблицы загружаются в виде GRanges аннотации в пределах указанных диапазонов GRanges.
library("rtracklayer")
 
# import  BED, GTF, WIG
setwd("~/Ubuntu One/evfr/redumeres/peak detection/")
bedpeak <- import("hs_merged_en10peaks.bed", asRangedData = F)
 
setwd("/media/data/db/danio/track/")
gtfanno <- import("zv9_danrer7_refgene.gtf", asRangedData = F)
 
setwd("~/Ubuntu One/evfr/redumeres/peak detection/")
wigcov <- import("dr1_uniq_peak-coverage.wig", asRangedData = F)
 
# UCSC browser is not cool. UCSC download tables is cool!
session <- browserSession("UCSC")
genome(session) <- "hg19"
track.names <- trackNames(ucscTableQuery(session)) # good list
# set track and table
loadingtrack <- "wgEncodeGencodeV17"
loadingtable <- "wgEncodeGencodePseudoGeneV17"
loadregion <- GRanges("chr10", IRanges(5000000, 8000000)
block <- getTable(
  ucscTableQuery(session, range = loadregion, track = loadingtrack, table = loadingtable)
)

Created by Pretty R at inside-R.org

R-Bioconductor - GenomicFeatures

Дисклеймер - немного по GenomicFeatures и TxDb



library("GenomicFeatures")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("mirbase.db")
library("FDb.UCSC.tRNAs")
 
# что вообще есть?
supportedUCSCtables() # много всего, но закачкой, сборкой и загрузкой заниматься не будем, все организмы слишком модельные
 
# базовые штуки
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene # переприсвоение для краткого имени
TranscriptsGR <- transcripts(txdb) # получение всех транскриптов в виде GRanges
ExonesGR <- exons(txdb) # экзоны
CDSGR <- cds(txdb, columns = c("tx_name", "cds_name")) # CDSы с доп. колонками
Prom <- promoters(txdb, upstream = 1000, downstream = 1000) # это не промотеры, а просто смещения
mirna <- microRNAs(txdb) # miRNA из MirBase
 
# выбор конкретных значений для колонок
vals <- list(tx_chrom = "chr5", tx_strand = "+") # значения ограничений
tx.sel.chr5 <- transcripts(txdb, vals = vals) # взята 5 хромосома и только + транскрипты
 
# группировка
GroupedTx <- transcriptsBy(txdb, by = "gene") # генерирует GRL по элементу на ген, в каждом элементе множество элементов гена
GroupedEx <- exonsBy(txdb, by = "tx") # GRL из траскриптов с экзонами
fiveUTR <- fiveUTRsByTranscript(txdb) #  типа 5'UTRы
intronsByTranscript(txdb) # типа наборы интронов
 
# выборка по пересечению
blocks <- GRanges(seqnames = rep("chr2", 10), # создаём пробный GRanges
  ranges = IRanges(start = c(1000, 5000, 15000, 20000, 30000, 50000, 100000, 120000, 130000, 200000),
    end = c(2000, 7000, 19000, 25000, 34000, 55000, 145000, 150000, 190000, 230000)),
  strand = rep("*", 10)
)
blockOverlap <- transcriptsByOverlaps(txdb, blocks) # выбираем что задевает его.
#### этого вполне хватит

Created by Pretty R at inside-R.org