Nothing
inst/misc/multiModels/aqp4-multi.R
In trena: Fit transcriptional regulatory networks using gene expression, priors, machine learning
library(TReNA)
library(RSQLite)
library(RPostgreSQL)
library(RUnit)
library(MotifDb); mdb <- MotifDb # convenience variable
library(splitstackshape)
#------------------------------------------------------------------------------------------------------------------------
fp.file <- "/Users/paul/s/work/priceLab/ohsu-aquaporin/aqp4.sqlite"
footprint.db.uri <- sprintf("sqlite://%s", fp.file);
fp.uri <- sprintf("sqlite://%s", fp.file)
genome.db.uri <- "postgres://bddsrds.globusgenomics.org/hg38"
#------------------------------------------------------------------------------------------------------------------------
tbl.mg <- read.table("~/github/TReNA/inst/extdata/motifGenes.tsv", sep="\t", as.is=TRUE, header=TRUE)
#------------------------------------------------------------------------------------------------------------------------
# if we need to do direct inspeciont
# db.fp <- dbConnect(dbDriver("SQLite"), fp.file)
#------------------------------------------------------------------------------------------------------------------------
tss <- 26865884 # minus strand
tssUpstream <- 500
tssDownstream <- 500
target.gene <- "AQP4"
#------------------------------------------------------------------------------------------------------------------------
if(!exists("mtx")){
load("~/github/projects/examples/microservices/trenaGeneModel/datasets/coryAD/rosmap_counts_matrix_normalized_geneSymbols_25031x638.RData")
#load(system.file(package="TReNA", "extdata/ampAD.154genes.mef2cTFs.278samples.RData"))
mtx <- asinh(mtx)
mtx.var <- apply(mtx, 1, var)
deleters <- which(mtx.var < 0.01)
if(length(deleters) > 0) # 15838 x 638
mtx <- mtx[-deleters,]
mtx <<- mtx
}
#------------------------------------------------------------------------------------------------------------------------
genome.name <- "hg38"
# rs3763043 disrupts MA0090.2
tbl.snp <- data.frame(target.gene=rep("AQP4", 5),
chromosome=rep("chr18", 5),
loc=c(26864410, 26865469, 26855623, 26855854, 26850565),
snp=c("rs3763040", "rs3875089", "rs335929", "rs3763043", "rs9951307"),
#shoulder=rep(2000, 5),
genome=rep(genome.name, 5),
stringsAsFactors=FALSE)
#------------------------------------------------------------------------------------------------------------------------
find.footprints <- function(chrom, start, end)
{
fp <- FootprintFinder(genome.db.uri, fp.uri, quiet=TRUE)
tbl.fp <- getFootprintsInRegion(fp, chrom, start, end)
} # find.footprints
#------------------------------------------------------------------------------------------------------------------------
find.snps.in.footprints <- function()
{
for(r in 1:nrow(tbl.snp)){
chrom <- tbl.snp$chromosome[r]
loc <- tbl.snp$loc[r]
shoulder <- 20
query <- sprintf("select * from regions where chrom='%s' and start < %d and endpos > %d",
chrom, loc+shoulder, loc-shoulder)
printf("---- %s: %d", tbl.snp$snp[r], tbl.snp$loc[r])
print(dbGetQuery(db.fp, query))
} # for r
} # find.snps.in.footprints
#------------------------------------------------------------------------------------------------------------------------
getFootprintCandidates <- function(target.gene, tssUpstream, tssDownstream)
{
regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- list()
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list();
filterSpec <- list(genomeDB=genome.db.uri,
footprintDB=footprint.db.uri,
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec)
filter <- FootprintFilter(filterSpec$genomeDB,
filterSpec$footprintDB,
filterSpec$geneCenteredSpec,
filterSpec$regionsSpec,
quiet=TRUE)
getCandidates(filter)
} # getFootprintCandidates
#------------------------------------------------------------------------------------------------------------------------
getDHSCandidates <- function(target.gene, tssUpstream, tssDownstream)
{
genome <- "hg38"
regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- NA_character_;
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list()
#variants <- NA_character_
variants <- "rs3875089"
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=75L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=85L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
quiet=FALSE))
getCandidates(filter)
} # getDHSCandidates
#------------------------------------------------------------------------------------------------------------------------
#createModel <- function(target.gene, tssUpstream, tssDownstream)
createModel <- function(target.gene, chrom, start, end, variants=NA_character_)
{
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list()
#regionsSpec <- list()
regionsSpec <- sprintf("%s:%d-%d", chrom, start, end)
filterSpec <- list(genomeDB=genome.db.uri,
footprintDB=footprint.db.uri,
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec)
filter <- FootprintFilter(filterSpec$genomeDB,
filterSpec$footprintDB,
filterSpec$geneCenteredSpec,
filterSpec$regionsSpec,
quiet=TRUE)
x.fp <- getCandidates(filter)
genome <- "hg38"
#regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- NA_character_;
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
#geneCenteredSpec <- list()
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=75L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=90L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
variants=variants,
quiet=FALSE))
x.dhs <- getCandidates(filter)
tfs.fp <- intersect(rownames(mtx), x.fp$tfs)
if(length(tfs.fp) > 0){
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.fp)
model.fp <- run(solver.wt)
}
else{
printf("no footprint-derived tfs")
x.fp <- list(tfs=NA_character_, tbl=data.frame())
model.fp <- list()
}
tfs.dhs <- intersect(rownames(mtx), x.dhs$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.dhs)
model.dhs <- run(solver.wt)
return(list(candidates.fp=x.fp, model.fp=model.fp, candidates.dhs=x.dhs, model.dhs=model.dhs))
} # createModel
#------------------------------------------------------------------------------------------------------------------------
createDHSModel <- function(target.gene, chrom, start, end, pwmMatchThreshold=80L, variants=NA_character_)
{
regionsSpec <- sprintf("%s:%d-%d", chrom, start, end)
geneCenteredSpec <- list()
genome <- "hg38"
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=pwmMatchThreshold,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=90L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
variants=variants,
quiet=FALSE))
x.dhs <- getCandidates(filter)
tfs.dhs <- intersect(rownames(mtx), x.dhs$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.dhs)
model.dhs <- run(solver.wt)
return(list(candidates=x.dhs, model=model.dhs))
} # createDHSModel
#------------------------------------------------------------------------------------------------------------------------
# 18:26865469 A/G
test.createModel <- function()
{
# chr18:26,860,742-26,882,685: includes all footprints and dhs clusters
# chr18:26,860,992-26,870,492: a 10kb region enriched with brain hint footprints
# chr18:26,864,303-26,866,143: a 1.8kb region including tss
start <- 26864303 #26860992 #26860742
end <- 26866143 #26870492 # 26882685
m1 <- createModel("AQP4", "chr18", start, end)
m1.mut <- createModel("AQP4", "chr18", start, end, variants="rs3875089")
# chr18:26865450-26865480: 30 bases around rs3875089
start <- 26865450
end <- 26865480
m2 <- createModel("AQP4", "chr18", start, end)
m2.mut <- createModel("AQP4", "chr18", start, end, variants="rs3875089")
if(length(m$model.fp) > 0){
head(m$model.dhs$edges)
# IncNodePurity gene.cor
# TEAD1 39.052237 0.7605383
# ATF7 19.679876 0.7163899
# SMAD9 15.450478 0.6694617
# SP3 13.706970 0.6837160
# NFE2L2 9.666957 0.6886920
# GLI2 9.514838 0.5750392
m$candidates.dhs$tbl[grep("TEAD", m$candidates.fp$tbl$tf),]
# chrom start end motifName length strand score1 score2 score3 tf
# 5851 chr18 26865890 26865897 MA0808.1 8 - 9 12.3265 6.51e-05 TEAD3;TEAD1;TEAD2;TEAD4
} # if fp results
head(m$model.dhs$edges)
m$candidates.dhs$tbl[grep("TEAD", m$candidates.dhs$tbl$tf),]
tbl.around.snp2 <- subset(m$candidates.dhs$tbl, motifStart<tbl.snp$loc[2] & motifEnd > tbl.snp$loc[2])
tfs.big <- rownames(subset(m$model.dhs$edges, IncNodePurity > 5))
# which of these high value tfs have a motif that includes snp2?
tbl.around.snp2[unlist(lapply(tfs.big, function(tf) grep(tf, tbl.around.snp2$tfs))),]
# motifName chrom motifStart motifEnd strand motifScore motifRelativeScore match regulatoryRegionStart regualtoryRegionEnd regulatorySequence variant tfs
# 17 MA0090.2 chr18 26865465 26865474 + 5.653203 0.8174939 AGCATCCCTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD1;TEAD2;TEAD3;TEAD4
# 111 MA0808.1 chr18 26865466 26865473 + 5.564402 0.8368466 GCATCCCT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD3;TEAD1;TEAD2;TEAD4
# 112 MA0809.1 chr18 26865465 26865474 + 5.470297 0.8060669 AGCATCCCTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD4;TEAD1;TEAD2;TEAD3
# 16 MA0084.1 chr18 26865468 26865476 - 4.535714 0.7016575 TTTGGCTAA 26865450 26865480 AGGATTTGGCTAAAAAG... wt SRY;SOX30;SOX15;SOX7;SOX17;SOX18;SOX8;SOX9;SOX10;SOX5;SOX6;SOX13;SOX4;SOX11;SOX12;SOX1;SOX2;SOX3;SOX14;SOX21
# 16.1 MA0084.1 chr18 26865468 26865476 - 4.535714 0.7016575 TTTGGCTAA 26865450 26865480 AGGATTTGGCTAAAAAG... wt SRY;SOX30;SOX15;SOX7;SOX17;SOX18;SOX8;SOX9;SOX10;SOX5;SOX6;SOX13;SOX4;SOX11;SOX12;SOX1;SOX2;SOX3;SOX14;SOX21
# 56 MA0597.1 chr18 26865467 26865475 + 4.070352 0.7155477 CATCCCTTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt THAP1;THAP10;THAP11;PRKRIR;THAP2;THAP3;THAP4;THAP5;THAP6;THAP7;THAP8;THAP9
# 130 MA0886.1 chr18 26865463 26865472 - 5.230989 0.7422321 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt EMX2;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 87 MA0710.1 chr18 26865463 26865472 - 5.096429 0.7505068 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt NOTO;NANOG;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 81 MA0699.1 chr18 26865463 26865472 - 4.596502 0.7295830 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt LBX2;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;BARHL1;BARHL2;NR2E1
# 127 MA0879.1 chr18 26865463 26865472 - 4.369845 0.7033855 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt DLX1;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 122 MA0876.1 chr18 26865464 26865471 - 4.064815 0.7189870 CTAAAAAG 26865450 26865480 AGGATTTGGCTAAAAAG... wt BSX;NANOG;NOTO;VENTX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# next up: do with the tbl.snp[2,] variant. does these tfs.big drop out?
} # test.createModel
#------------------------------------------------------------------------------------------------------------------------
createModels <- function(spec, mtx)
{
stopifnot(all(c("target.gene", "chromosome", "loc", "snp", "shoulder", "genome") %in% names(spec)))
regionsSpec <- with(spec, sprintf("%s:%d-%d", chromosome, loc-shoulder, loc+shoulder))
dhsFilterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=spec$genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=80L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=list(),
regionsSpec=c(regionsSpec),
variants=spec$snp)
hdf.wt <- with(dhsFilterSpec,
HumanDHSFilter(genomeName=genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=pwmMatchPercentageThreshold,
geneCenteredSpec=geneCenteredSpec,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
quiet=TRUE))
x.wt <- getCandidates(hdf.wt)
# find any motif matches, in regulatory regions, which contain the variant site
subset(x.wt$tbl, motifStart < spec$loc & motifEnd > spec$loc)
hdf.mut <- with(dhsFilterSpec,
HumanDHSFilter(genomeName=genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=pwmMatchPercentageThreshold,
geneCenteredSpec=geneCenteredSpec,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
variants=variants,
quiet=TRUE))
x.mut <- getCandidates(hdf.mut)
# find any motif matches, in regulatory regions, which contain the variant site
#subset(x.wt$tbl, motifStart <= spec$loc & motifEnd >= spec$loc)
#subset(x.mut$tbl, motifStart <= spec$loc & motifEnd >= spec$loc)
tfs.wt <- intersect(rownames(mtx), x.wt$tfs)
tfs.mut <- intersect(rownames(mtx), x.mut$tfs)
identical.tfs <- (length(tfs.wt) == length(tfs.mut)) & all(tfs.wt %in% tfs.mut)
printf("%d wt tfs, %d mut tfs, identical? %s", length(tfs.wt), length(tfs.mut), identical.tfs)
model.wt <- data.frame()
model.mut <- data.frame()
if(length(tfs.wt) > 0){
solver.wt <- RandomForestSolver(mtx, targetGene=spec$target.gene, candidateRegulators=tfs.wt)
model.wt <- run(solver.wt)
}
else{
printf("no candidate tfs for wt")
}
if(length(tfs.mut) > 0){
solver.mut <- RandomForestSolver(mtx, targetGene=spec$target.gene, candidateRegulators=tfs.mut)
model.mut <- run(solver.mut)
}
else{
printf("no candidate tfs for wt")
}
list(wt.candidates=x.wt, mut.candidates=x.mut, wt.model=model.wt, mut.model=model.mut)
} # createModels
#------------------------------------------------------------------------------------------------------------------------
explore.aqp4 <- function()
{
start <- 26864303 #26860992 #26860742
end <- 26866143 #26870492 # 26882685
start.300bp <- 26865352
end.300bp <- 26865619
start <- start.300bp
end <- end.300bp
#m1 <- createDHSModel("AQP4", "chr18", start, end, 95L)
tv <- TReNA.Viz(portRange=11011:11051)
# addBedTrackFromHostedFile(tv,
# trackName="brain HINT",
# uri="http://pshannon.systemsbiology.net/annotations/brain_hint.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/brain_hint.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
# addBedTrackFromHostedFile(tv,
# trackName="skin HINT",
# uri="http://pshannon.systemsbiology.net/annotations/skin_hint.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/skin_hint.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
addBedTrackFromHostedFile(tv,
trackName="EncodeDHSclustered",
uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz",
index.uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz.tbi",
displayMode="SQUISHED",
color="blue")
addBedTrackFromHostedFile(tv,
trackName="AQP4 snps",
uri="http://pshannon.systemsbiology.net/annotations/aqp4-chr18-snps.bed",
index.uri=NA,
displayMode="SQUISHED",
color="red")
#tbl.bed <- m1$candidates$tbl[, c("chrom", "motifStart", "motifEnd", "motifName", "motifScore")]
#addBedTrackFromDataFrame(tv, "DHS motifs", tbl.bed, color="green")
tv
} # explore.aqp4
#----------------------------------------------------------------------------------------------------
displayAllTracks <- function(tv, chrom, start, end)
{
if(!exists("db"))
db <<- dbConnect(PostgreSQL(), user= "trena", password="trena", host="whovian")
dbNames <- c("brain_hint", "brain_wellington",
"lymphoblast_hint", "lymphoblast_wellington",
"skin_hint", "skin_wellington")
actual.dbNames <- grep("_", dbGetQuery(db, "select datname from pg_database")[, 1], value=TRUE)
dbNames <- intersect(dbNames, actual.dbNames)
dbNames <- list()
genome.db.uri <- "postgres://whovian/hg38"
# chrom <- "chr18"
# # based on all 5 snps
# start.17kb.all.snps <- 26849000
# end.17kb.all.snps <- 26866000
# # based on tss
# start <- tss - 2000
# end <- tss + 2000
# start <- 26847227
# end <- 26869162
# start.81kb <- 26824211
# end.81kb <- 26905511
# start <- start.81kb
# end <- end.81kb
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
tbl.regions <- data.frame();
load("aqp4_fp.RData")
tbl.cory <- subset(fp$AQP4$tbl, start >= start & endpos <= end)[, c(3, 4, 5, 1, 15)]
tbl.cory$source <- "cory"
tbl.regions <- rbind(tbl.regions, tbl.cory)
addBedTrackFromDataFrame(tv, "coryFP", tbl.cory, displayMode="COLLAPSED", color="cyan")
for(dbName in dbNames){
printf("--- querying %s", dbName)
fp.db.uri <- sprintf("postgres://whovian/%s", dbName)
fpf <- FootprintFinder(genome.db.uri, fp.db.uri)
tbl.fp <- getFootprintsInRegion(fpf, chrom, start, end)
printf("%40s: %5d footprints", dbName, nrow(tbl.fp))
if(nrow(tbl.fp) > 0){
tbl.bed <- tbl.fp[, c("chrom", "start", "endpos", "name", "score2")]
tbl.bed$source <- dbName
tbl.regions <- rbind(tbl.regions, tbl.bed)
addBedTrackFromDataFrame(tv, dbName, tbl.bed, displayMode="COLLAPSED", color="green")
}
closeDatabaseConnections(fpf)
} # for dbName
addBedTrackFromDataFrame(tv, "fp.all", tbl.regions, color="orange")
save(tbl.regions, file="all.aqp4.footprints.RData")
invisible(tbl.regions)
} # displayAllTracks
#----------------------------------------------------------------------------------------------------
getCandidates <- function(tbl.regions)
{
i <- grep("endpos", colnames(tbl.regions))
if(length(i) > 0)
colnames(tbl.regions)[i] <- "end"
motifMatcher <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
#x <- findMatchesByChromosomalRegion(motifMatcher, tbl.regions, pwmMatchMinimumAsPercentage=92)
tbl.regions.single <- data.frame(chrom="chr18", start=26863884, end=26867884, stringsAsFactors=FALSE)
x92 = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=92)
x70 = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=70)
print(system.time(x2 <- findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=72)))
tfs <- intersect(rownames(mtx), x92$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
model.wt <- run(solver.wt)
snps <- tbl.snp$snp
x92.snps = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=92, variants=snps)
} # getCandidates
#----------------------------------------------------------------------------------------------------
newFriday.rs3875089 <- function()
{
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
rsid <- "rs3875089"
loc <- subset(tbl.snp, snp==rsid)$loc
tbl.regions.noSeq <- data.frame(chrom="chr18",
start=loc - 10,
end=loc + 10,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
tbl.mut <- getSequence(mm, tbl.regions.noSeq, rsid)
checkTrue(tbl.wt$seq != tbl.mut$seq)
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80, rsid)
} # newFriday.rs3875089
#----------------------------------------------------------------------------------------------------
identifyPerturbedMotifs <- function(chrom="chr18", start=26850565, end=26865469)
{
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=min(tbl.snp$loc) - 100,
end=max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
tbl.mut <- getSequence(mm, tbl.regions.noSeq, tbl.snp$snp)
checkEquals(nchar(tbl.wt$seq), nchar(tbl.mut$seq))
for(i in 1:nchar(tbl.mut$seq)){
wt.base <- substr(tbl.wt$seq, i, i)
mut.base <- substr(tbl.mut$seq, i, i)
if(wt.base != mut.base)
printf("variant at position %d: %s -> %s", i, wt.base, mut.base)
} # for i
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=70)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=70, variants=tbl.snp$snp)
tbl.wt.freq <- as.data.frame(table(x.wt$tbl$motifName))
tbl.wt.freq <- tbl.wt.freq[order(tbl.wt.freq$Freq, decreasing=TRUE),]
colnames(tbl.wt.freq) <- c("motif", "wtCount")
tbl.mut.freq <- as.data.frame(table(x.mut$tbl$motifName))
tbl.mut.freq <- tbl.mut.freq[order(tbl.mut.freq$Freq, decreasing=TRUE),]
colnames(tbl.mut.freq) <- c("motif", "mutCount")
# "MA0751.1" lost in
tbl.counts <- merge(tbl.wt.freq, tbl.mut.freq, by="motif", all=TRUE);
tbl.counts$mutCount[is.na(tbl.counts$mutCount)] <- 0
tbl.counts$wtCount[is.na(tbl.counts$wtCount)] <- 0
tbl.counts$diff <- tbl.counts$wtCount - tbl.counts$mutCount
tbl.counts <- tbl.counts[order(abs(tbl.counts$diff), decreasing=TRUE),]
tbl.mg <- read.table("~/github/TReNA/inst/extdata/motifGenes.tsv", sep="\t", as.is=TRUE, header=TRUE)
tbl.counts$motif <- as.character(tbl.counts$motif)
tfs <- unlist(lapply(tbl.counts$motif, function(m) paste(subset(tbl.mg, motif==m)$tf.gene, collapse=";")))
tbl.counts$tfs <- tfs
mapMotifToTF <- function(motif){
motif.stem <- strsplit(motif, ".", fixed=TRUE)[[1]][1]
motif.stem.regex <- sprintf("%s\\.", motif.stem)
tfs <- unique(mcols(query(mdb, motif.stem.regex))$geneSymbol)
paste(tfs, collapse=";")
}
mdb.tfs <- unlist(lapply((tbl.counts$motif), mapMotifToTF))
tbl.counts$motifDB.tfs <- mdb.tfs
save(tbl.counts, file="aqp4.tbl.counts.5kbUpstream.RData")
wt.ma0090.locs <- unique(with(subset(x.wt$tbl, motifName=="MA0090.2"), sprintf("%s:%d-%d", chrom, motifStart, motifEnd)))
mut.ma0090.locs <- unique(with(subset(x.mut$tbl, motifName=="MA0090.2"), sprintf("%s:%d-%d", chrom, motifStart, motifEnd)))
setdiff(wt.ma0090.locs, mut.ma0090.locs) # [1] "chr18:26855847-26855856" "chr18:26865465-26865474"
subset(x.wt$tbl, (motifStart==26855847 | motifStart==26865465) & motifName=="MA0090.2")
subset(x.mut$tbl, (motifStart==26855847 | motifStart==26865465) & motifName=="MA0090.2")
# motifName chrom motifStart motifEnd strand motifScore motifRelativeScore match chromStart chromEnd seq status tf
# |
# 3462 MA0090.2 chr18 26855847 26855856 + 5.734231 0.8292111 AATATTCCAG 26850465 26865569 CCCATATATATGCTCACAATTGATAATTATTCTAATG... wt TEAD1;TEAD2;TEAD3;TEAD4
# |
# 3481 MA0090.2 chr18 26865465 26865474 + 5.653203 0.8174939 AGCATCCCTT 26850465 26865569 CCCATATATATGCTCACAATTGATAATTATTCTAATG... wt TEAD1;TEAD2;TEAD3;TEAD4
# | |
# A [ 218 734 212 1132 0 265 29 0 621 280 ]
# C [ 454 100 920 63 10 33 1132 1132 109 346 ]
# G [ 147 398 21 0 4 0 2 1 206 138 ]
# T [ 314 51 18 32 1132 867 24 291 511 367 ]
tbl.snp
# target.gene chromosome loc snp shoulder genome
# AQP4 chr18 26864410 rs3763040 2000 hg38
# AQP4 chr18 26865469 rs3875089 2000 hg38 A->G (reverse)
# AQP4 chr18 26855623 rs335929 2000 hg38
# AQP4 chr18 26855854 rs3763043 2000 hg38 A->G (reverse?)
# AQP4 chr18 26850565 rs9951307 2000 hg38
} # identifyPerturbedMotifs
#----------------------------------------------------------------------------------------------------
buildModel <- function()
{
if(!exists("tbl.count"))
load("aqp4.tbl.counts.5kbUpstream.RData", envir=.GlobalEnv)
printf("motifs mapped to genes in motifDB: %d", length(intersect(toupper(tbl.counts$motifDB.tfs), rownames(mtx)))) # 101
tfs.oldStyle <- unique(unlist(strsplit((tbl.counts$tfs), ";")))
printf("tfs.oldStyle: %d", length(tfs.oldStyle)) # 714
printf("motifs mapped to genes old style: %d", length(intersect(tfs.oldStyle, rownames(mtx)))) # 466
tfs <- intersect(rownames(mtx), tfs.oldStyle)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
} # buildModel
#----------------------------------------------------------------------------------------------------
init.tv <- function(all.tracks=FALSE)
{
tv <- TReNA.Viz(portRange=11011:11051)
#
# addBedTrackFromHostedFile(tv,
# trackName="EncodeDHSclustered",
# uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
#
#
# addBedTrackFromHostedFile(tv,
# trackName="AQP4 snps",
# uri="http://pshannon.systemsbiology.net/annotations/aqp4-chr18-snps.bed",
# index.uri=NA,
# displayMode="SQUISHED",
# color="red")
# showGenomicRegion(tv, "chr18:26,865,340-26,865,658")
# if(all.tracks){
# displayAllTracks()
# } # all.tracks
tv
} # init.tv
#----------------------------------------------------------------------------------------------------
runBasic <- function(chrom="chr18", start=26865462, end=26865867, min.motif.score=80L)
{
chrom="chr18"; start=26865462; end=26865867; min.motif.score=80L # short and simple, just 400bp
chrom="chr18"; start=26849820; end=26871166 # 21kb
tv <- init.tv()
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=min.motif.score)
tfs.oldStyle <- unique(unlist(strsplit((x.wt$tbl$tf), ";")))
printf("tfs.oldStyle: %d", length(tfs.oldStyle)) # 714
printf("motifs mapped to genes old style: %d", length(intersect(tfs.oldStyle, rownames(mtx)))) # 466
tfs <- intersect(rownames(mtx), tfs.oldStyle)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
#addBedTrackFromDataFrame(tv, "all motifs", x.wt$tbl[, c(2,3,4,1, 6)], color="grey", displayMode="COLLAPSED")
for(tf.strong in tfs.strong){
trackName <- sprintf("%s motifs", tf.strong)
#addBedTrackFromDataFrame(tv, trackName, x.wt$tbl[grep(tf.strong, x.wt$tbl$tf), c(2,3,4,1, 6)], color="green")
}
# expand x.wt$tbl (motifs & their tfs) to one row per motif/tf. then remove all but tfs.strong
# the resulting table should be ready for TReNA-Viz::geneRegulatoryModelToGraph
printf("----- expanding tfs")
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
printf("converting model to graph, %d tfs, %d motif sites", nrow(tbl.model), nrow(tbl.motifs))
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
xyz <- 99
} # runBasic
#----------------------------------------------------------------------------------------------------
runWithAugmentedFootprints <- function()
{
start.81kb <- 26824211
end.81kb <- 26905511
chrom <- "chr18"
start.300bp <- 26865352
end.300bp <- 26865619
start.1425bp <- 26864212
end.1425bp <- 26865636
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tbl.regions <- displayAllTracks(chrom, start, end)
dim(tbl.regions) # 88938 x 5
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=83)
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
tv <- init.tv()
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
#--- now remove TEAD1 motifs, rebuild model
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants="rs3875089")
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.mut$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.mut <- run(solver.mut)
tfs.strong <- rownames(subset(model.mut$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.mut$edges)
tbl.model <- data.frame(tf=rownames(model.mut$edges),
randomForest=model.mut$edges$IncNodePurity,
pearson=model.mut$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
tbl.model2 <- subset(tbl.model, tf != "TEAD1")
tbl.motifs2 <- subset(tbl.motifs, tf != "TEAD1")
system.time(g2 <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model2, tbl.motifs2))
printf("adding model layout");
system.time(g2.lo <- TReNA:::addGeneModelLayout(g2))
printf("addGraph")
addGraph(tv, g2.lo)
loadStyle(tv, "style.js")
save(tbl.model, tbl.motifs, tbl.model2, tbl.motifs2, g.lo, g2.lo, file="aqp4.300bp.modelAndGraph.RData")
} # runWithAugmentedFootprints
#----------------------------------------------------------------------------------------------------
# AQP4.fp.upstream.00207.L10.MA0090.2: 26865677 26865686 - 5.651008 0.8171764 TTGGAGTGTT
# AQP4.fp.upstream.00199.L10.MA0090.2: 26865685 26865694 + 5.562714 0.8044084 CTCATGCCTT
quick <- function()
{
chrom <- "chr18"
start.300bp <- 26865450
end.300bp <- 26865480
start.tss.m300 <- tss - 300
end.tss <- tss
start <- start.tss.m300
end <- end.tss
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tv <- init.tv(all.tracks=FALSE)
# tbl.regions <- displayAllTracks(tv, chrom, start, end)
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
dim(tbl.regions) # 88938 x 5
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80)
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
chrom <- "chr18"
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tv2 <- init.tv(all.tracks=FALSE)
#tbl.regions <- displayAllTracks(tv2, chrom, start, end)
showGenomicRegion(tv2, sprintf("%s:%d-%d", chrom, start, end))
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants="rs3875089")
tbl.trimmed <- subset(x.mut$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.mut)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv2, g.lo)
loadStyle(tv2, "style.js")
#tbl.model2 <- subset(tbl.model, tf != "TEAD1")
#tbl.motifs2 <- subset(tbl.motifs, tf != "TEAD1")
system.time(g2 <- geneRegulatoryModelToGraph(tv2, "AQP4", tbl.model2, tbl.motifs2))
printf("adding model layout");
system.time(g2.lo <- TReNA:::addGeneModelLayout(g2))
printf("addGraph")
addGraph(tv2, g2.lo)
loadStyle(tv2, "style.js")
xyz <- 99
} # quick
#------------------------------------------------------------------------------------------------------------------------
buildModelForRegion <- function(chrom, start, end, variants=NA)
{
#tv <- init.tv(all.tracks=FALSE)
#showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
x <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants=variants)
tbl.trimmed <- subset(x$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.regulatoryRegions <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.regulatoryRegions))
tbl.regulatoryRegions$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.regulatoryRegions$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.mut)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.regulatoryRegions <- subset(tbl.regulatoryRegions, tf %in% tfs.strong)
tbl.regulatoryRegions$distance.from.tss <- tss - tbl.regulatoryRegions$motifStart
regulatoryRegion.names <- unlist(lapply(1:nrow(tbl.regulatoryRegions), function(i){
distance.from.tss <- tbl.regulatoryRegions$distance.from.tss[i]
region.size <- nchar(tbl.regulatoryRegions$match[i])
motif.name <- tbl.regulatoryRegions$motifName[i]
if(distance.from.tss < 0)
sprintf("%s.fp.downstream.%05d.L%d.%s", target.gene, abs(distance.from.tss), region.size, motif.name)
else
sprintf("%s.fp.upstream.%05d.L%d.%s", target.gene, abs(distance.from.tss), region.size, motif.name)
}))
tbl.regulatoryRegions$regionName <- regulatoryRegion.names
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
list(tbl.model=tbl.model, tbl.regulatoryRegions=tbl.regulatoryRegions)
} # buildModelForRegion
#------------------------------------------------------------------------------------------------------------------------
demo <- function()
{
pad<-8
i<-2
x.wt <- buildModelForRegion("chr18", tbl.snp$loc[i] - pad, tbl.snp$loc[i] + pad)
x.mut <- buildModelForRegion("chr18", tbl.snp$loc[i] - pad, tbl.snp$loc[i] + pad, tbl.snp$snp[2])
save(x.wt, x.mut, file="../../extdata/twoAQP4modelsForTesting.RData")
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.regulatoryRegions))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
} # demo
#------------------------------------------------------------------------------------------------------------------------
# wt <-
tead1.motifs <- unique(subset(tbl.mg, tf.gene=="TEAD1")$motif)
Try the trena package in your browser
Any scripts or data that you put into this service are public.
trena documentation built on Nov. 15, 2020, 2:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(TReNA)
library(RSQLite)
library(RPostgreSQL)
library(RUnit)
library(MotifDb); mdb <- MotifDb # convenience variable
library(splitstackshape)
#------------------------------------------------------------------------------------------------------------------------
fp.file <- "/Users/paul/s/work/priceLab/ohsu-aquaporin/aqp4.sqlite"
footprint.db.uri <- sprintf("sqlite://%s", fp.file);
fp.uri <- sprintf("sqlite://%s", fp.file)
genome.db.uri <- "postgres://bddsrds.globusgenomics.org/hg38"
#------------------------------------------------------------------------------------------------------------------------
tbl.mg <- read.table("~/github/TReNA/inst/extdata/motifGenes.tsv", sep="\t", as.is=TRUE, header=TRUE)
#------------------------------------------------------------------------------------------------------------------------
# if we need to do direct inspeciont
# db.fp <- dbConnect(dbDriver("SQLite"), fp.file)
#------------------------------------------------------------------------------------------------------------------------
tss <- 26865884 # minus strand
tssUpstream <- 500
tssDownstream <- 500
target.gene <- "AQP4"
#------------------------------------------------------------------------------------------------------------------------
if(!exists("mtx")){
load("~/github/projects/examples/microservices/trenaGeneModel/datasets/coryAD/rosmap_counts_matrix_normalized_geneSymbols_25031x638.RData")
#load(system.file(package="TReNA", "extdata/ampAD.154genes.mef2cTFs.278samples.RData"))
mtx <- asinh(mtx)
mtx.var <- apply(mtx, 1, var)
deleters <- which(mtx.var < 0.01)
if(length(deleters) > 0) # 15838 x 638
mtx <- mtx[-deleters,]
mtx <<- mtx
}
#------------------------------------------------------------------------------------------------------------------------
genome.name <- "hg38"
# rs3763043 disrupts MA0090.2
tbl.snp <- data.frame(target.gene=rep("AQP4", 5),
chromosome=rep("chr18", 5),
loc=c(26864410, 26865469, 26855623, 26855854, 26850565),
snp=c("rs3763040", "rs3875089", "rs335929", "rs3763043", "rs9951307"),
#shoulder=rep(2000, 5),
genome=rep(genome.name, 5),
stringsAsFactors=FALSE)
#------------------------------------------------------------------------------------------------------------------------
find.footprints <- function(chrom, start, end)
{
fp <- FootprintFinder(genome.db.uri, fp.uri, quiet=TRUE)
tbl.fp <- getFootprintsInRegion(fp, chrom, start, end)
} # find.footprints
#------------------------------------------------------------------------------------------------------------------------
find.snps.in.footprints <- function()
{
for(r in 1:nrow(tbl.snp)){
chrom <- tbl.snp$chromosome[r]
loc <- tbl.snp$loc[r]
shoulder <- 20
query <- sprintf("select * from regions where chrom='%s' and start < %d and endpos > %d",
chrom, loc+shoulder, loc-shoulder)
printf("---- %s: %d", tbl.snp$snp[r], tbl.snp$loc[r])
print(dbGetQuery(db.fp, query))
} # for r
} # find.snps.in.footprints
#------------------------------------------------------------------------------------------------------------------------
getFootprintCandidates <- function(target.gene, tssUpstream, tssDownstream)
{
regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- list()
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list();
filterSpec <- list(genomeDB=genome.db.uri,
footprintDB=footprint.db.uri,
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec)
filter <- FootprintFilter(filterSpec$genomeDB,
filterSpec$footprintDB,
filterSpec$geneCenteredSpec,
filterSpec$regionsSpec,
quiet=TRUE)
getCandidates(filter)
} # getFootprintCandidates
#------------------------------------------------------------------------------------------------------------------------
getDHSCandidates <- function(target.gene, tssUpstream, tssDownstream)
{
genome <- "hg38"
regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- NA_character_;
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list()
#variants <- NA_character_
variants <- "rs3875089"
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=75L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=85L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
quiet=FALSE))
getCandidates(filter)
} # getDHSCandidates
#------------------------------------------------------------------------------------------------------------------------
#createModel <- function(target.gene, tssUpstream, tssDownstream)
createModel <- function(target.gene, chrom, start, end, variants=NA_character_)
{
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
geneCenteredSpec <- list()
#regionsSpec <- list()
regionsSpec <- sprintf("%s:%d-%d", chrom, start, end)
filterSpec <- list(genomeDB=genome.db.uri,
footprintDB=footprint.db.uri,
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec)
filter <- FootprintFilter(filterSpec$genomeDB,
filterSpec$footprintDB,
filterSpec$geneCenteredSpec,
filterSpec$regionsSpec,
quiet=TRUE)
x.fp <- getCandidates(filter)
genome <- "hg38"
#regionsSpec <- "chr18:26865459-26865479"
#regionsSpec <- NA_character_;
#geneCenteredSpec <- list(targetGene=target.gene, tssUpstream=tssUpstream, tssDownstream=tssDownstream)
#geneCenteredSpec <- list()
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=75L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=90L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
variants=variants,
quiet=FALSE))
x.dhs <- getCandidates(filter)
tfs.fp <- intersect(rownames(mtx), x.fp$tfs)
if(length(tfs.fp) > 0){
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.fp)
model.fp <- run(solver.wt)
}
else{
printf("no footprint-derived tfs")
x.fp <- list(tfs=NA_character_, tbl=data.frame())
model.fp <- list()
}
tfs.dhs <- intersect(rownames(mtx), x.dhs$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.dhs)
model.dhs <- run(solver.wt)
return(list(candidates.fp=x.fp, model.fp=model.fp, candidates.dhs=x.dhs, model.dhs=model.dhs))
} # createModel
#------------------------------------------------------------------------------------------------------------------------
createDHSModel <- function(target.gene, chrom, start, end, pwmMatchThreshold=80L, variants=NA_character_)
{
regionsSpec <- sprintf("%s:%d-%d", chrom, start, end)
geneCenteredSpec <- list()
genome <- "hg38"
filterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=pwmMatchThreshold,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=geneCenteredSpec,
regionsSpec=regionsSpec,
variants=variants)
filter <- with(filterSpec, HumanDHSFilter(genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=90L,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
geneCenteredSpec=geneCenteredSpec,
variants=variants,
quiet=FALSE))
x.dhs <- getCandidates(filter)
tfs.dhs <- intersect(rownames(mtx), x.dhs$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene=target.gene, candidateRegulators=tfs.dhs)
model.dhs <- run(solver.wt)
return(list(candidates=x.dhs, model=model.dhs))
} # createDHSModel
#------------------------------------------------------------------------------------------------------------------------
# 18:26865469 A/G
test.createModel <- function()
{
# chr18:26,860,742-26,882,685: includes all footprints and dhs clusters
# chr18:26,860,992-26,870,492: a 10kb region enriched with brain hint footprints
# chr18:26,864,303-26,866,143: a 1.8kb region including tss
start <- 26864303 #26860992 #26860742
end <- 26866143 #26870492 # 26882685
m1 <- createModel("AQP4", "chr18", start, end)
m1.mut <- createModel("AQP4", "chr18", start, end, variants="rs3875089")
# chr18:26865450-26865480: 30 bases around rs3875089
start <- 26865450
end <- 26865480
m2 <- createModel("AQP4", "chr18", start, end)
m2.mut <- createModel("AQP4", "chr18", start, end, variants="rs3875089")
if(length(m$model.fp) > 0){
head(m$model.dhs$edges)
# IncNodePurity gene.cor
# TEAD1 39.052237 0.7605383
# ATF7 19.679876 0.7163899
# SMAD9 15.450478 0.6694617
# SP3 13.706970 0.6837160
# NFE2L2 9.666957 0.6886920
# GLI2 9.514838 0.5750392
m$candidates.dhs$tbl[grep("TEAD", m$candidates.fp$tbl$tf),]
# chrom start end motifName length strand score1 score2 score3 tf
# 5851 chr18 26865890 26865897 MA0808.1 8 - 9 12.3265 6.51e-05 TEAD3;TEAD1;TEAD2;TEAD4
} # if fp results
head(m$model.dhs$edges)
m$candidates.dhs$tbl[grep("TEAD", m$candidates.dhs$tbl$tf),]
tbl.around.snp2 <- subset(m$candidates.dhs$tbl, motifStart<tbl.snp$loc[2] & motifEnd > tbl.snp$loc[2])
tfs.big <- rownames(subset(m$model.dhs$edges, IncNodePurity > 5))
# which of these high value tfs have a motif that includes snp2?
tbl.around.snp2[unlist(lapply(tfs.big, function(tf) grep(tf, tbl.around.snp2$tfs))),]
# motifName chrom motifStart motifEnd strand motifScore motifRelativeScore match regulatoryRegionStart regualtoryRegionEnd regulatorySequence variant tfs
# 17 MA0090.2 chr18 26865465 26865474 + 5.653203 0.8174939 AGCATCCCTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD1;TEAD2;TEAD3;TEAD4
# 111 MA0808.1 chr18 26865466 26865473 + 5.564402 0.8368466 GCATCCCT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD3;TEAD1;TEAD2;TEAD4
# 112 MA0809.1 chr18 26865465 26865474 + 5.470297 0.8060669 AGCATCCCTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt TEAD4;TEAD1;TEAD2;TEAD3
# 16 MA0084.1 chr18 26865468 26865476 - 4.535714 0.7016575 TTTGGCTAA 26865450 26865480 AGGATTTGGCTAAAAAG... wt SRY;SOX30;SOX15;SOX7;SOX17;SOX18;SOX8;SOX9;SOX10;SOX5;SOX6;SOX13;SOX4;SOX11;SOX12;SOX1;SOX2;SOX3;SOX14;SOX21
# 16.1 MA0084.1 chr18 26865468 26865476 - 4.535714 0.7016575 TTTGGCTAA 26865450 26865480 AGGATTTGGCTAAAAAG... wt SRY;SOX30;SOX15;SOX7;SOX17;SOX18;SOX8;SOX9;SOX10;SOX5;SOX6;SOX13;SOX4;SOX11;SOX12;SOX1;SOX2;SOX3;SOX14;SOX21
# 56 MA0597.1 chr18 26865467 26865475 + 4.070352 0.7155477 CATCCCTTT 26865450 26865480 AGGATTTGGCTAAAAAG... wt THAP1;THAP10;THAP11;PRKRIR;THAP2;THAP3;THAP4;THAP5;THAP6;THAP7;THAP8;THAP9
# 130 MA0886.1 chr18 26865463 26865472 - 5.230989 0.7422321 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt EMX2;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 87 MA0710.1 chr18 26865463 26865472 - 5.096429 0.7505068 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt NOTO;NANOG;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 81 MA0699.1 chr18 26865463 26865472 - 4.596502 0.7295830 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt LBX2;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;BARHL1;BARHL2;NR2E1
# 127 MA0879.1 chr18 26865463 26865472 - 4.369845 0.7033855 GCTAAAAAGC 26865450 26865480 AGGATTTGGCTAAAAAG... wt DLX1;NANOG;NOTO;VENTX;BSX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# 122 MA0876.1 chr18 26865464 26865471 - 4.064815 0.7189870 CTAAAAAG 26865450 26865480 AGGATTTGGCTAAAAAG... wt BSX;NANOG;NOTO;VENTX;HHEX;HLX;EN1;EN2;EMX1;EMX2;DLX1;DLX2;DLX3;DLX4;DLX5;DLX6;DBX1;DBX2;VAX1;VAX2;TLX1;TLX2;TLX3;BARX1;BARX2;HMX1;HMX2;HMX3;MSX1;MSX2;LBX1;LBX2;BARHL1;BARHL2;NR2E1
# next up: do with the tbl.snp[2,] variant. does these tfs.big drop out?
} # test.createModel
#------------------------------------------------------------------------------------------------------------------------
createModels <- function(spec, mtx)
{
stopifnot(all(c("target.gene", "chromosome", "loc", "snp", "shoulder", "genome") %in% names(spec)))
regionsSpec <- with(spec, sprintf("%s:%d-%d", chromosome, loc-shoulder, loc+shoulder))
dhsFilterSpec <- list(filterType="EncodeDNaseClusters",
genomeName=spec$genome,
encodeTableName="wgEncodeRegDnaseClustered",
pwmMatchPercentageThreshold=80L,
geneInfoDB="postgres://whovian/gtf",
geneCenteredSpec=list(),
regionsSpec=c(regionsSpec),
variants=spec$snp)
hdf.wt <- with(dhsFilterSpec,
HumanDHSFilter(genomeName=genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=pwmMatchPercentageThreshold,
geneCenteredSpec=geneCenteredSpec,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
quiet=TRUE))
x.wt <- getCandidates(hdf.wt)
# find any motif matches, in regulatory regions, which contain the variant site
subset(x.wt$tbl, motifStart < spec$loc & motifEnd > spec$loc)
hdf.mut <- with(dhsFilterSpec,
HumanDHSFilter(genomeName=genomeName,
encodeTableName=encodeTableName,
pwmMatchPercentageThreshold=pwmMatchPercentageThreshold,
geneCenteredSpec=geneCenteredSpec,
geneInfoDatabase.uri=geneInfoDB,
regionsSpec=regionsSpec,
variants=variants,
quiet=TRUE))
x.mut <- getCandidates(hdf.mut)
# find any motif matches, in regulatory regions, which contain the variant site
#subset(x.wt$tbl, motifStart <= spec$loc & motifEnd >= spec$loc)
#subset(x.mut$tbl, motifStart <= spec$loc & motifEnd >= spec$loc)
tfs.wt <- intersect(rownames(mtx), x.wt$tfs)
tfs.mut <- intersect(rownames(mtx), x.mut$tfs)
identical.tfs <- (length(tfs.wt) == length(tfs.mut)) & all(tfs.wt %in% tfs.mut)
printf("%d wt tfs, %d mut tfs, identical? %s", length(tfs.wt), length(tfs.mut), identical.tfs)
model.wt <- data.frame()
model.mut <- data.frame()
if(length(tfs.wt) > 0){
solver.wt <- RandomForestSolver(mtx, targetGene=spec$target.gene, candidateRegulators=tfs.wt)
model.wt <- run(solver.wt)
}
else{
printf("no candidate tfs for wt")
}
if(length(tfs.mut) > 0){
solver.mut <- RandomForestSolver(mtx, targetGene=spec$target.gene, candidateRegulators=tfs.mut)
model.mut <- run(solver.mut)
}
else{
printf("no candidate tfs for wt")
}
list(wt.candidates=x.wt, mut.candidates=x.mut, wt.model=model.wt, mut.model=model.mut)
} # createModels
#------------------------------------------------------------------------------------------------------------------------
explore.aqp4 <- function()
{
start <- 26864303 #26860992 #26860742
end <- 26866143 #26870492 # 26882685
start.300bp <- 26865352
end.300bp <- 26865619
start <- start.300bp
end <- end.300bp
#m1 <- createDHSModel("AQP4", "chr18", start, end, 95L)
tv <- TReNA.Viz(portRange=11011:11051)
# addBedTrackFromHostedFile(tv,
# trackName="brain HINT",
# uri="http://pshannon.systemsbiology.net/annotations/brain_hint.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/brain_hint.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
# addBedTrackFromHostedFile(tv,
# trackName="skin HINT",
# uri="http://pshannon.systemsbiology.net/annotations/skin_hint.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/skin_hint.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
addBedTrackFromHostedFile(tv,
trackName="EncodeDHSclustered",
uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz",
index.uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz.tbi",
displayMode="SQUISHED",
color="blue")
addBedTrackFromHostedFile(tv,
trackName="AQP4 snps",
uri="http://pshannon.systemsbiology.net/annotations/aqp4-chr18-snps.bed",
index.uri=NA,
displayMode="SQUISHED",
color="red")
#tbl.bed <- m1$candidates$tbl[, c("chrom", "motifStart", "motifEnd", "motifName", "motifScore")]
#addBedTrackFromDataFrame(tv, "DHS motifs", tbl.bed, color="green")
tv
} # explore.aqp4
#----------------------------------------------------------------------------------------------------
displayAllTracks <- function(tv, chrom, start, end)
{
if(!exists("db"))
db <<- dbConnect(PostgreSQL(), user= "trena", password="trena", host="whovian")
dbNames <- c("brain_hint", "brain_wellington",
"lymphoblast_hint", "lymphoblast_wellington",
"skin_hint", "skin_wellington")
actual.dbNames <- grep("_", dbGetQuery(db, "select datname from pg_database")[, 1], value=TRUE)
dbNames <- intersect(dbNames, actual.dbNames)
dbNames <- list()
genome.db.uri <- "postgres://whovian/hg38"
# chrom <- "chr18"
# # based on all 5 snps
# start.17kb.all.snps <- 26849000
# end.17kb.all.snps <- 26866000
# # based on tss
# start <- tss - 2000
# end <- tss + 2000
# start <- 26847227
# end <- 26869162
# start.81kb <- 26824211
# end.81kb <- 26905511
# start <- start.81kb
# end <- end.81kb
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
tbl.regions <- data.frame();
load("aqp4_fp.RData")
tbl.cory <- subset(fp$AQP4$tbl, start >= start & endpos <= end)[, c(3, 4, 5, 1, 15)]
tbl.cory$source <- "cory"
tbl.regions <- rbind(tbl.regions, tbl.cory)
addBedTrackFromDataFrame(tv, "coryFP", tbl.cory, displayMode="COLLAPSED", color="cyan")
for(dbName in dbNames){
printf("--- querying %s", dbName)
fp.db.uri <- sprintf("postgres://whovian/%s", dbName)
fpf <- FootprintFinder(genome.db.uri, fp.db.uri)
tbl.fp <- getFootprintsInRegion(fpf, chrom, start, end)
printf("%40s: %5d footprints", dbName, nrow(tbl.fp))
if(nrow(tbl.fp) > 0){
tbl.bed <- tbl.fp[, c("chrom", "start", "endpos", "name", "score2")]
tbl.bed$source <- dbName
tbl.regions <- rbind(tbl.regions, tbl.bed)
addBedTrackFromDataFrame(tv, dbName, tbl.bed, displayMode="COLLAPSED", color="green")
}
closeDatabaseConnections(fpf)
} # for dbName
addBedTrackFromDataFrame(tv, "fp.all", tbl.regions, color="orange")
save(tbl.regions, file="all.aqp4.footprints.RData")
invisible(tbl.regions)
} # displayAllTracks
#----------------------------------------------------------------------------------------------------
getCandidates <- function(tbl.regions)
{
i <- grep("endpos", colnames(tbl.regions))
if(length(i) > 0)
colnames(tbl.regions)[i] <- "end"
motifMatcher <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
#x <- findMatchesByChromosomalRegion(motifMatcher, tbl.regions, pwmMatchMinimumAsPercentage=92)
tbl.regions.single <- data.frame(chrom="chr18", start=26863884, end=26867884, stringsAsFactors=FALSE)
x92 = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=92)
x70 = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=70)
print(system.time(x2 <- findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=72)))
tfs <- intersect(rownames(mtx), x92$tfs)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
model.wt <- run(solver.wt)
snps <- tbl.snp$snp
x92.snps = findMatchesByChromosomalRegion(motifMatcher, tbl.regions.single, pwmMatchMinimumAsPercentage=92, variants=snps)
} # getCandidates
#----------------------------------------------------------------------------------------------------
newFriday.rs3875089 <- function()
{
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
rsid <- "rs3875089"
loc <- subset(tbl.snp, snp==rsid)$loc
tbl.regions.noSeq <- data.frame(chrom="chr18",
start=loc - 10,
end=loc + 10,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
tbl.mut <- getSequence(mm, tbl.regions.noSeq, rsid)
checkTrue(tbl.wt$seq != tbl.mut$seq)
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80, rsid)
} # newFriday.rs3875089
#----------------------------------------------------------------------------------------------------
identifyPerturbedMotifs <- function(chrom="chr18", start=26850565, end=26865469)
{
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=FALSE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=min(tbl.snp$loc) - 100,
end=max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
tbl.mut <- getSequence(mm, tbl.regions.noSeq, tbl.snp$snp)
checkEquals(nchar(tbl.wt$seq), nchar(tbl.mut$seq))
for(i in 1:nchar(tbl.mut$seq)){
wt.base <- substr(tbl.wt$seq, i, i)
mut.base <- substr(tbl.mut$seq, i, i)
if(wt.base != mut.base)
printf("variant at position %d: %s -> %s", i, wt.base, mut.base)
} # for i
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=70)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=70, variants=tbl.snp$snp)
tbl.wt.freq <- as.data.frame(table(x.wt$tbl$motifName))
tbl.wt.freq <- tbl.wt.freq[order(tbl.wt.freq$Freq, decreasing=TRUE),]
colnames(tbl.wt.freq) <- c("motif", "wtCount")
tbl.mut.freq <- as.data.frame(table(x.mut$tbl$motifName))
tbl.mut.freq <- tbl.mut.freq[order(tbl.mut.freq$Freq, decreasing=TRUE),]
colnames(tbl.mut.freq) <- c("motif", "mutCount")
# "MA0751.1" lost in
tbl.counts <- merge(tbl.wt.freq, tbl.mut.freq, by="motif", all=TRUE);
tbl.counts$mutCount[is.na(tbl.counts$mutCount)] <- 0
tbl.counts$wtCount[is.na(tbl.counts$wtCount)] <- 0
tbl.counts$diff <- tbl.counts$wtCount - tbl.counts$mutCount
tbl.counts <- tbl.counts[order(abs(tbl.counts$diff), decreasing=TRUE),]
tbl.mg <- read.table("~/github/TReNA/inst/extdata/motifGenes.tsv", sep="\t", as.is=TRUE, header=TRUE)
tbl.counts$motif <- as.character(tbl.counts$motif)
tfs <- unlist(lapply(tbl.counts$motif, function(m) paste(subset(tbl.mg, motif==m)$tf.gene, collapse=";")))
tbl.counts$tfs <- tfs
mapMotifToTF <- function(motif){
motif.stem <- strsplit(motif, ".", fixed=TRUE)[[1]][1]
motif.stem.regex <- sprintf("%s\\.", motif.stem)
tfs <- unique(mcols(query(mdb, motif.stem.regex))$geneSymbol)
paste(tfs, collapse=";")
}
mdb.tfs <- unlist(lapply((tbl.counts$motif), mapMotifToTF))
tbl.counts$motifDB.tfs <- mdb.tfs
save(tbl.counts, file="aqp4.tbl.counts.5kbUpstream.RData")
wt.ma0090.locs <- unique(with(subset(x.wt$tbl, motifName=="MA0090.2"), sprintf("%s:%d-%d", chrom, motifStart, motifEnd)))
mut.ma0090.locs <- unique(with(subset(x.mut$tbl, motifName=="MA0090.2"), sprintf("%s:%d-%d", chrom, motifStart, motifEnd)))
setdiff(wt.ma0090.locs, mut.ma0090.locs) # [1] "chr18:26855847-26855856" "chr18:26865465-26865474"
subset(x.wt$tbl, (motifStart==26855847 | motifStart==26865465) & motifName=="MA0090.2")
subset(x.mut$tbl, (motifStart==26855847 | motifStart==26865465) & motifName=="MA0090.2")
# motifName chrom motifStart motifEnd strand motifScore motifRelativeScore match chromStart chromEnd seq status tf
# |
# 3462 MA0090.2 chr18 26855847 26855856 + 5.734231 0.8292111 AATATTCCAG 26850465 26865569 CCCATATATATGCTCACAATTGATAATTATTCTAATG... wt TEAD1;TEAD2;TEAD3;TEAD4
# |
# 3481 MA0090.2 chr18 26865465 26865474 + 5.653203 0.8174939 AGCATCCCTT 26850465 26865569 CCCATATATATGCTCACAATTGATAATTATTCTAATG... wt TEAD1;TEAD2;TEAD3;TEAD4
# | |
# A [ 218 734 212 1132 0 265 29 0 621 280 ]
# C [ 454 100 920 63 10 33 1132 1132 109 346 ]
# G [ 147 398 21 0 4 0 2 1 206 138 ]
# T [ 314 51 18 32 1132 867 24 291 511 367 ]
tbl.snp
# target.gene chromosome loc snp shoulder genome
# AQP4 chr18 26864410 rs3763040 2000 hg38
# AQP4 chr18 26865469 rs3875089 2000 hg38 A->G (reverse)
# AQP4 chr18 26855623 rs335929 2000 hg38
# AQP4 chr18 26855854 rs3763043 2000 hg38 A->G (reverse?)
# AQP4 chr18 26850565 rs9951307 2000 hg38
} # identifyPerturbedMotifs
#----------------------------------------------------------------------------------------------------
buildModel <- function()
{
if(!exists("tbl.count"))
load("aqp4.tbl.counts.5kbUpstream.RData", envir=.GlobalEnv)
printf("motifs mapped to genes in motifDB: %d", length(intersect(toupper(tbl.counts$motifDB.tfs), rownames(mtx)))) # 101
tfs.oldStyle <- unique(unlist(strsplit((tbl.counts$tfs), ";")))
printf("tfs.oldStyle: %d", length(tfs.oldStyle)) # 714
printf("motifs mapped to genes old style: %d", length(intersect(tfs.oldStyle, rownames(mtx)))) # 466
tfs <- intersect(rownames(mtx), tfs.oldStyle)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
} # buildModel
#----------------------------------------------------------------------------------------------------
init.tv <- function(all.tracks=FALSE)
{
tv <- TReNA.Viz(portRange=11011:11051)
#
# addBedTrackFromHostedFile(tv,
# trackName="EncodeDHSclustered",
# uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz",
# index.uri="http://pshannon.systemsbiology.net/annotations/dhsClusters_hg38.bed.gz.tbi",
# displayMode="SQUISHED",
# color="blue")
#
#
# addBedTrackFromHostedFile(tv,
# trackName="AQP4 snps",
# uri="http://pshannon.systemsbiology.net/annotations/aqp4-chr18-snps.bed",
# index.uri=NA,
# displayMode="SQUISHED",
# color="red")
# showGenomicRegion(tv, "chr18:26,865,340-26,865,658")
# if(all.tracks){
# displayAllTracks()
# } # all.tracks
tv
} # init.tv
#----------------------------------------------------------------------------------------------------
runBasic <- function(chrom="chr18", start=26865462, end=26865867, min.motif.score=80L)
{
chrom="chr18"; start=26865462; end=26865867; min.motif.score=80L # short and simple, just 400bp
chrom="chr18"; start=26849820; end=26871166 # 21kb
tv <- init.tv()
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=min.motif.score)
tfs.oldStyle <- unique(unlist(strsplit((x.wt$tbl$tf), ";")))
printf("tfs.oldStyle: %d", length(tfs.oldStyle)) # 714
printf("motifs mapped to genes old style: %d", length(intersect(tfs.oldStyle, rownames(mtx)))) # 466
tfs <- intersect(rownames(mtx), tfs.oldStyle)
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
#addBedTrackFromDataFrame(tv, "all motifs", x.wt$tbl[, c(2,3,4,1, 6)], color="grey", displayMode="COLLAPSED")
for(tf.strong in tfs.strong){
trackName <- sprintf("%s motifs", tf.strong)
#addBedTrackFromDataFrame(tv, trackName, x.wt$tbl[grep(tf.strong, x.wt$tbl$tf), c(2,3,4,1, 6)], color="green")
}
# expand x.wt$tbl (motifs & their tfs) to one row per motif/tf. then remove all but tfs.strong
# the resulting table should be ready for TReNA-Viz::geneRegulatoryModelToGraph
printf("----- expanding tfs")
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
printf("converting model to graph, %d tfs, %d motif sites", nrow(tbl.model), nrow(tbl.motifs))
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
xyz <- 99
} # runBasic
#----------------------------------------------------------------------------------------------------
runWithAugmentedFootprints <- function()
{
start.81kb <- 26824211
end.81kb <- 26905511
chrom <- "chr18"
start.300bp <- 26865352
end.300bp <- 26865619
start.1425bp <- 26864212
end.1425bp <- 26865636
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tbl.regions <- displayAllTracks(chrom, start, end)
dim(tbl.regions) # 88938 x 5
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=83)
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
tv <- init.tv()
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
#--- now remove TEAD1 motifs, rebuild model
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants="rs3875089")
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.mut$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.mut <- run(solver.mut)
tfs.strong <- rownames(subset(model.mut$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.mut$edges)
tbl.model <- data.frame(tf=rownames(model.mut$edges),
randomForest=model.mut$edges$IncNodePurity,
pearson=model.mut$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
tbl.model2 <- subset(tbl.model, tf != "TEAD1")
tbl.motifs2 <- subset(tbl.motifs, tf != "TEAD1")
system.time(g2 <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model2, tbl.motifs2))
printf("adding model layout");
system.time(g2.lo <- TReNA:::addGeneModelLayout(g2))
printf("addGraph")
addGraph(tv, g2.lo)
loadStyle(tv, "style.js")
save(tbl.model, tbl.motifs, tbl.model2, tbl.motifs2, g.lo, g2.lo, file="aqp4.300bp.modelAndGraph.RData")
} # runWithAugmentedFootprints
#----------------------------------------------------------------------------------------------------
# AQP4.fp.upstream.00207.L10.MA0090.2: 26865677 26865686 - 5.651008 0.8171764 TTGGAGTGTT
# AQP4.fp.upstream.00199.L10.MA0090.2: 26865685 26865694 + 5.562714 0.8044084 CTCATGCCTT
quick <- function()
{
chrom <- "chr18"
start.300bp <- 26865450
end.300bp <- 26865480
start.tss.m300 <- tss - 300
end.tss <- tss
start <- start.tss.m300
end <- end.tss
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tv <- init.tv(all.tracks=FALSE)
# tbl.regions <- displayAllTracks(tv, chrom, start, end)
showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
dim(tbl.regions) # 88938 x 5
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
tbl.wt <- getSequence(mm, tbl.regions.noSeq)
printf("findMatchesByChromosomalRegion, size: %d", 1 + end - start);
x.wt <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=80)
#x.fp <- findMatchesByChromosomalRegion(mm, tbl.regions, pwmMatchMinimumAsPercentage=75)
tbl.trimmed <- subset(x.wt$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.wt <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.wt)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
chrom <- "chr18"
start <- tbl.snp$loc[2] - 10
end <- tbl.snp$loc[2] + 10
tv2 <- init.tv(all.tracks=FALSE)
#tbl.regions <- displayAllTracks(tv2, chrom, start, end)
showGenomicRegion(tv2, sprintf("%s:%d-%d", chrom, start, end))
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
x.mut <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants="rs3875089")
tbl.trimmed <- subset(x.mut$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.motifs <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.motifs))
tbl.motifs$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.motifs$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.mut)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.motifs <- subset(tbl.motifs, tf %in% tfs.strong)
tbl.motifs$distance.from.tss <- tss - tbl.motifs$motifStart
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.motifs))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv2, g.lo)
loadStyle(tv2, "style.js")
#tbl.model2 <- subset(tbl.model, tf != "TEAD1")
#tbl.motifs2 <- subset(tbl.motifs, tf != "TEAD1")
system.time(g2 <- geneRegulatoryModelToGraph(tv2, "AQP4", tbl.model2, tbl.motifs2))
printf("adding model layout");
system.time(g2.lo <- TReNA:::addGeneModelLayout(g2))
printf("addGraph")
addGraph(tv2, g2.lo)
loadStyle(tv2, "style.js")
xyz <- 99
} # quick
#------------------------------------------------------------------------------------------------------------------------
buildModelForRegion <- function(chrom, start, end, variants=NA)
{
#tv <- init.tv(all.tracks=FALSE)
#showGenomicRegion(tv, sprintf("%s:%d-%d", chrom, start, end))
mm <- MotifMatcher(name="mm", genomeName="hg38", quiet=TRUE)
tbl.regions.noSeq <- data.frame(chrom=chrom,
start=start, # min(tbl.snp$loc) - 100,
end=end, #max(tbl.snp$loc) + 5000,
stringsAsFactors=FALSE)
x <- findMatchesByChromosomalRegion(mm, tbl.regions.noSeq, pwmMatchMinimumAsPercentage=75, variants=variants)
tbl.trimmed <- subset(x$tbl, nchar(tf) != 0)
tfs.split <- strsplit(tbl.trimmed$tf, ";")
length(tfs.split) # [1] 36929
counts <- unlist(lapply(tfs.split, length))
tfs.split.vec <- unlist(tfs.split)
tbl.regulatoryRegions <- expandRows(tbl.trimmed, counts, count.is.col=FALSE, drop=FALSE)
checkEquals(length(tfs.split.vec), nrow(tbl.regulatoryRegions))
tbl.regulatoryRegions$tf <- tfs.split.vec
tfs <- intersect(rownames(mtx), unique(tbl.regulatoryRegions$tf))
solver.mut <- RandomForestSolver(mtx, targetGene="AQP4", candidateRegulators=tfs)
printf("calling randomForestSolver");
model.wt <- run(solver.mut)
tfs.strong <- rownames(subset(model.wt$edges, IncNodePurity > 2))
tbl.regulatoryRegions <- subset(tbl.regulatoryRegions, tf %in% tfs.strong)
tbl.regulatoryRegions$distance.from.tss <- tss - tbl.regulatoryRegions$motifStart
regulatoryRegion.names <- unlist(lapply(1:nrow(tbl.regulatoryRegions), function(i){
distance.from.tss <- tbl.regulatoryRegions$distance.from.tss[i]
region.size <- nchar(tbl.regulatoryRegions$match[i])
motif.name <- tbl.regulatoryRegions$motifName[i]
if(distance.from.tss < 0)
sprintf("%s.fp.downstream.%05d.L%d.%s", target.gene, abs(distance.from.tss), region.size, motif.name)
else
sprintf("%s.fp.upstream.%05d.L%d.%s", target.gene, abs(distance.from.tss), region.size, motif.name)
}))
tbl.regulatoryRegions$regionName <- regulatoryRegion.names
printf("----- expanding tfs, done")
count <- nrow(model.wt$edges)
tbl.model <- data.frame(tf=rownames(model.wt$edges),
randomForest=model.wt$edges$IncNodePurity,
pearson=model.wt$edges$gene.cor,
spearman=rep(0, count),
betaLasso=rep(0, count),
pcaMax=rep(0, count),
concordance=rep(0, count),
stringsAsFactors=FALSE)
tbl.model <- subset(tbl.model, randomForest >= 2)
list(tbl.model=tbl.model, tbl.regulatoryRegions=tbl.regulatoryRegions)
} # buildModelForRegion
#------------------------------------------------------------------------------------------------------------------------
demo <- function()
{
pad<-8
i<-2
x.wt <- buildModelForRegion("chr18", tbl.snp$loc[i] - pad, tbl.snp$loc[i] + pad)
x.mut <- buildModelForRegion("chr18", tbl.snp$loc[i] - pad, tbl.snp$loc[i] + pad, tbl.snp$snp[2])
save(x.wt, x.mut, file="../../extdata/twoAQP4modelsForTesting.RData")
system.time(g <- geneRegulatoryModelToGraph(tv, "AQP4", tbl.model, tbl.regulatoryRegions))
printf("adding model layout");
system.time(g.lo <- TReNA:::addGeneModelLayout(g))
printf("addGraph")
addGraph(tv, g.lo)
loadStyle(tv, "style.js")
} # demo
#------------------------------------------------------------------------------------------------------------------------
# wt <-
tead1.motifs <- unique(subset(tbl.mg, tf.gene=="TEAD1")$motif)
Try the trena package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.