inst/doc/vtpnet.R

In vtpnet: variant-transcription factor-phenotype networks

### R code from vignette source 'vtpnet.Rnw'

###################################################
### code chunk number 1: lkco
###################################################
#
# tags formatted for display
#
diseaseTags = c("Ankylosing\\\nspondylitis", "Asthma",
      "Celiac\\\ndisease", "Crohn's\\\ndisease",
      "Multiple\\\nsclerosis", "Primary\\\nbiliary\\\ncirrhosis",
      "Psoriasis", "Rheumatoid\\\narthritis", 
      "Systemic\\\nlupus\\\nerythematosus",
      "Systemic\\\nsclerosis", "Type 1\\\ndiabetes", 
      "Ulcerative\\\ncolitis"
)
TFtags = c("ELF3", "MEF2A", "TCF3", "PAX4", "STAT3",
   "ESR1", "POU2F1", "STAT1", "YY1", "SP1", "CDC5L",
   "NR3C1", "EGR1", "PPARG", "HNF4A", "REST", "PPARA",
   "AR", "NFKB1", "HNF1A", "TFAP2A")
# define adjacency matrix
adjm = matrix(1, nr=length(diseaseTags), nc=length(TFtags))
dimnames(adjm) = list(diseaseTags, TFtags)
library(graph)
cvtp = ugraph(aM2bpG(adjm)) # complete (V)TP network; variants not involved yet


###################################################
### code chunk number 2: dog1
###################################################
library(Rgraphviz)
#flat = function(x, g) c(x, edges(g)[[x]])
#sub = subGraph(unique(c(flat("Crohn's\\\ndisease", cvtp), 
#    flat("Ulcerative\\\ncolitis", cvtp))), cvtp)
sub = subGraph(unique(c(diseaseTags[1:4], TFtags[1:6])), cvtp)
plot(sub, attrs=list(node=list(shape="box", fixedsize=FALSE)))
#plot(cvtp, attrs=list(graph=list(margin=c(.5,.5), size=c(4.1,4.1)), 
#   node=list(shape="box", fixedsize=FALSE, height=1)))


###################################################
### code chunk number 3: lkvd
###################################################
library(vtpnet)
data(maurGWAS)
length(maurGWAS)
names(values(maurGWAS))


###################################################
### code chunk number 4: lkpax
###################################################
data(pax4)
length(pax4)
pax4[1:4]


###################################################
### code chunk number 5: lksearch (eval = FALSE)
###################################################
## library(foreach)
## library(doParallel)
## registerDoParallel(cores=12)
## library(BSgenome.Hsapiens.UCSC.hg19)
## library(MotifDb)
## sn = seqnames(Hsapiens)[1:24]
## pax4 = query(MotifDb, "pax4")[[1]]
## ans = foreach(i=1:24) %dopar% {
##  cat(i)
##  subj = Hsapiens[[ sn[i] ]]
##  matchPWM( pax4, subj, "75%" )
## }
## pax4_75 = 
##  do.call(c, lapply(1:length(ans), function(x) 
##   {GRanges(sn[x], as(ans[[x]], "IRanges"))}))
## save(pax4_75, file="pax4_75.rda")


###################################################
### code chunk number 6: lkmat
###################################################
data(pax4_85)
vp_pax4_85 = maurGWAS[ overlapsAny(maurGWAS, pax4_85) ]
length(vp_pax4_85)
data(pax4_75)
vp_pax4_75 = maurGWAS[ overlapsAny(maurGWAS, pax4_75) ]
length(vp_pax4_75)


###################################################
### code chunk number 7: lkfim
###################################################
vp_pax4_fimo = maurGWAS[ overlapsAny(maurGWAS, pax4) ]
length(vp_pax4_fimo)


###################################################
### code chunk number 8: lkta
###################################################
u75 = unique(vp_pax4_75$disease_trait)
ufimo = unique(vp_pax4_fimo$disease_trait)
length(setdiff(u75, ufimo))
length(setdiff(ufimo, u75))


###################################################
### code chunk number 9: domap
###################################################
data(cancerMap)
requireNamespace("gwascat")
load(system.file("legacy/gwrngs19.rda", package="gwascat"))
cangw = filterGWASbyMap( gwrngs19, cancerMap )
getOneHits( pax4, cangw, "fimo" )


###################################################
### code chunk number 10: algiz
###################################################
altize = function(chtag = "21",
#
# from sketch by Herve Pages, May 2013
#
  slpack="SNPlocs.Hsapiens.dbSNP.20120608",
  hgpack ="BSgenome.Hsapiens.UCSC.hg19",
  faElFun = function(x) sub("%%TAG%%", x, "alt%%TAG%%chr"),
  faTargFun = function(x) 
     sub("%%TAG%%", x, "alt%%TAG%%_hg19.fa")) {
    require(slpack, character.only=TRUE)
    require(hgpack, character.only=TRUE)
    require("ShortRead", character.only=TRUE)
    chk = grep("ch|chr", chtag)
    if (length(chk)>0) {
      warning("clearing prefix ch or chr from chtag")
      chtag = gsub("ch|chr", "", chtag)
      }
    snpgettag = paste0("ch", chtag)
    ggettag = paste0("chr", chtag)
    cursnps = getSNPlocs(snpgettag, as.GRanges=TRUE)
    curgenome = unmasked(Hsapiens[[ggettag]])
    ref_allele = 
     strsplit(as.character(curgenome[start(cursnps)]),
        NULL, fixed=TRUE)[[1L]]
    all_alleles = IUPAC_CODE_MAP[cursnps$alleles_as_ambig]
    alt_alleles = mapply( function(ref,all)
      sub(ref, "", all, fixed=TRUE),
        ref_allele, all_alleles, USE.NAMES=FALSE)
    cursnps$ref_allele = ref_allele
    cursnps$alt_alleles = alt_alleles
    cursnps$one_alt = substr(cursnps$alt_alleles, 1, 1)
    altg = list(replaceLetterAt(curgenome, start(cursnps),
      cursnps$one_alt))
    names(altg) = faElFun(chtag)
    writeFasta(DNAStringSet(altg), file=faTargFun(chtag))
}


###################################################
### code chunk number 11: lksess
###################################################
sessionInfo()