katwre

> head(cage)
GRanges object with 6 ranges and 1 metadata column:
      seqnames               ranges strand |       tpm
         <Rle>            <IRanges>  <Rle> | <numeric>
  [1]    chr22 [42017258, 42017373]      + |  534.6991
  [2]    chr22 [24236548, 24236627]      + |  640.9699
  [3]    chr21 [34915313, 34915452]      + |  461.0562
  [4]    chr22 [18111391, 18111584]      - |  403.6615
  [5]    chr22 [38071637, 38071688]      + |  450.2464
  [6]    chr22 [36783979, 36784055]      - |  373.9197

katwre

# ---------------------------------------------------------------------------- #
#' Get combination of elements of a character vector
#' 
#' @vec a vector of characters
find.combinations = function(vec){
  
  all.comb = expand.grid(vec,
                         vec,
                         KEEP.OUT.ATTRS = FALSE,
                         stringsAsFactors = FALSE)

katwre

getTabixByChr2gr = function(tbxFile, chr, return.type){
  require(Rsamtools)
  require(methylKit)
  require(GenomicRanges)
  from = methylKit:::getTabixByChr(tbxFile,chr=chr,
                        return.type=return.type)
  GRanges(seqnames=as.character(from$V1),
        ranges=IRanges(start=from$V2, end=from$V3),
        strand=from$V4, from[,5:ncol(from)])  
}

katwre

## this function tries to fetch the reference genes for the given assembly
fetchTable <- function(table.name,
                       table.loc = NULL,
                       assembly) {
  
  if(is.null(table.loc)) table.loc <- paste0(table.name, assembly,".bed")
  
  ## import local bed file if available
  if( file.exists(table.loc) ) {
    ## parse it 

katwre

# http://www.vincebuffalo.com/blog/2012/03/12/using-bioconductor-to-analyze-your-23andme-data.html

library(gwascat)
data(ebicat37) 
library(GenomicRanges)

data="~/projects/23andme/genome_v4.txt"

katwre

uniform_tfbs = data.frame(cell.line=c("GM12878" ,"H1-hESC" ,"K562", "HeLa-S3", "HepG2", "HUVEC", "A549", "IMR-90",	"MCF-7"),
                          schema=c("wgEncodeAwgTfbsUtaGm12878CtcfUniPk", "wgEncodeAwgTfbsUtaH1hescCtcfUniPk", "wgEncodeAwgTfbsUtaK562CtcfUniPk", "wgEncodeAwgTfbsUtaHelas3CtcfUniPk",
                                   "wgEncodeAwgTfbsUtaHepg2CtcfUniPk", "wgEncodeAwgTfbsUtaHuvecCtcfUniPk", "wgEncodeAwgTfbsHaibA549Ctcfsc5916Pcr1xDex100nmUniPk",
                                   "wgEncodeAwgTfbsSydhImr90CtcfbIggrabUniPk", "wgEncodeAwgTfbsUtaMcf7CtcfSerumstimUniPk"),
                          stringsAsFactors = FALSE)


#' UCSC table to GRanges object (only chrom, start and end)
#' 
#' @param table UCSC table

katwre

library(RcppRoll)
library(microbenchmark)

x <- runif(1000,0.5,1.5)
microbenchmark(
   runmed(x = x, k = 3),
    roll_median(x, 3),
   times=1000L
)

katwre

bampath = "example.bam"
command=paste0("samtools idxstats ", bampath," | awk 'BEGIN {a=0} {a += $3 } END{print a }'")
tot.mapped=as.numeric(try(system(command,intern=TRUE)))

katwre

destdir = "~/" #directory where fastq files will be stored
nmbr_of_cores = 20
SRAdb_path = '~/SRAmetadb.sqlite'
exp_ids=c("ERX620541")

require(SRAdb)
sqlfile <- SRAdb_path
if(!file.exists(SRAdb_path)) sqlfile <- getSRAdbFile()
sra_con <- dbConnect(SQLite(),sqlfile) #conection to SRA database

katwre

   target.paired.end = GRanges(rep(c(1),each=7), 
                              IRanges(c(7,7,8,9,9,24,25), width=39),
                              strand=rep(c("-"), times=7))
  windows.paired.end = GRanges(rep(c(1),each=3), IRanges(c(7,8,24), width=10), 
                               strand=rep("-", times=3))

   plotir <- function(ir,i, col) { arrows(start(ir)-.5,i,end(ir)+.5,i,code=3,angle=90,lwd=3, col=col) }
  plot(0,0,xlim=c(0,65),ylim=c(0,11),type="n",xlab="",ylab="",xaxt="n")
  axis(1,0:65)
  abline(v=0:65 + .5,col=rgb(0,0,0,.5))