## ----style, echo = FALSE, results = 'asis'-------------------------------------------------------- BiocStyle::markdown() options(width=100, max.print=1000) knitr::opts_chunk$set( eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")), cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")), error=FALSE) ## ----txdb----------------------------------------------------------------------------------------- suppressPackageStartupMessages({ library("TxDb.Rnorvegicus.UCSC.rn5.refGene") }) txdb <- TxDb.Rnorvegicus.UCSC.rn5.refGene ## find all genes ratGenes <- genes(txdb) ## list all sequences seqinfo(ratGenes) ## subset to contain only standard chromosomes ratGenes <- keepStandardChromosomes(ratGenes) ## find gene 'Acsl5' acsl5 <- ratGenes[which(mcols(ratGenes)$gene_id==94340),] acsl5 ## ----bsgenome------------------------------------------------------------------------------------- suppressPackageStartupMessages({ library(BSgenome.Rnorvegicus.UCSC.rn5) }) ratSeq <- BSgenome.Rnorvegicus.UCSC.rn5 class(ratSeq) ## get the sequence acsl5_sequence <- getSeq(ratSeq, acsl5) ## calculate the GC content letterFrequency(acsl5_sequence, "GC", as.prob=TRUE) ## ----select-rat----------------------------------------------------------------------------------- library("Rattus.norvegicus") ## get a mapping between all entrex id and gene names ratGeneNames <- select(Rattus.norvegicus, ratGenes$gene_id, columns=c("SYMBOL", 'GENEID'), keytype="GENEID") ## match the entrz id with result before subsetting idx <- match(ratGeneNames$GENEID, ratGenes$gene_id) ## add mactched result to GRanges mcols(ratGenes) <- ratGeneNames[idx,] ratGenes ## ----gtf-to-gr, eval=FALSE------------------------------------------------------------------------ # # library(AnnotationHub) # ah = AnnotationHub() # # ## find the file # gtf_humans <- query(ah , c("gtf", "Homo sapiens", "grch38","80")) # gtf_humans # # ## download the file # gtfFile <- ah[["AH47066"]] # # ## create a txdb # library(GenomicFeatures) # txdb <- makeTxDbFromGRanges(gtfFile) #may take some time.. # txdb # # ## get the genes from the taxdb object # humanGenes <- genes(txdb) ## ----chainfile------------------------------------------------------------------------------------ ## a) get the chain file ## load the package and query the files to find the file we want library(AnnotationHub) ah = AnnotationHub() query(ah , c("rattus", "rn5", "rn6")) ## learn more about the file you want ah["AH14761"] ## download the file ratChain <- ah[["AH14761"]] ratChain ## b) perform the liftover library(rtracklayer) lft <- liftOver(acsl5, ratChain) lft ## ----sessionInfo---------------------------------------------------------------------------------- sessionInfo()