library(RUnit)
library(snpFoot);
library(PrivateCoryData)
library(igvR)
#------------------------------------------------------------------------------------------------------------------------
Sys.setlocale("LC_ALL", "C")
#------------------------------------------------------------------------------------------------------------------------
printf <- function(...) print(noquote(sprintf(...)))
#------------------------------------------------------------------------------------------------------------------------
if(!exists("tbl.fpAnnotated")){
printf("loading tbl.fpAnnotated from PrivateCoryData package...")
data(tbl.fpAnnotated)
}
if(!exists("tbl.fpAnnotated")){
printf("loading tbl.fpAnnotated from PrivateCoryData package...")
data(tbl.fpAnnotated)
}
if(!exists("tbl.motifToMultipleGenes")){
printf("loading tbl.motifToMultipleGenes from PrivateCoryData package...")
data(tbl.motifToMultipleGenes)
}
#------------------------------------------------------------------------------------------------------------------------
runTests <- function()
{
test_1_footprint_SAT2B_tf()
test_1_footprint_two_tfs()
test_emptyResults()
test_360k_chr5_alzheimers_region()
test_findSNPsInFootprints()
test_findSNPsNearFootprints()
test_LXH1_matching()
test_tfGrabber()
test_tfGrabber.new()
test_run.tfGrabber.new()
test_intersectingLocs()
test_displayGWAS()
test_intersectMarietSnpsWithGWAS()
test_displayAllEncodeFootprints()
test_runLift.hg19to38()
test_selectBestAmongDuplicates()
} # runTests
#------------------------------------------------------------------------------------------------------------------------
# tbl.motifToMultipleGenes (from seth's "motif_to_tf_mappings_with_tfclass_include_multiple.csv" tells us that
# that SATB2 binds to pwm motif MA0679.1.
# strategy:
# 1) find a footprint in tbl.fpAnnotated by direct inspection: grep("MA0679", tbl.fpAnnotated$info)[1] # 891
# 2) use the chromosomal location of that cherry-picked footprint
# 3) make sure that findTFsInFootprints retrieves that footprint, and annotates it to SATB2
test_1_footprint_SAT2B_tf <- function()
{
printf("--- test_1_footprint_SAT2B_tf")
loc <- tbl.fpAnnotated[grep("MA0679", tbl.fpAnnotated$motifName)[1], c("chr", "mfpStart", "mfpEnd")]
target.tf <- "SATB2"
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$mfpStart, endLoc=loc$mfpEnd,
transcriptionFactors=target.tf)
# tbl.fpAnnotated identifies 3 motifs in this one footprint, and each of them is associated with SATB2
# we therefore expect 3 rows, each describing the same footprint, and each with a different motif
checkEquals(dim(tbl), c(3, 18))
checkEquals(tbl$tfsMatched, rep("SATB2", 3))
checkEquals(tbl$motifName, c("MA0679.1", "MA0756.1", "MA0757.1"))
} # test_1_footprint_SAT2B_tf
#------------------------------------------------------------------------------------------------------------------------
test_1_footprint_two_tfs <- function()
{
printf("--- test_1_footprint_two_tfs")
loc <- tbl.fpAnnotated[grep("MA0679", tbl.fpAnnotated$motifName)[1], c("chr", "mfpStart", "mfpEnd")]
target.tf <- c("SATB2", "CUX1")
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$mfpStart, endLoc=loc$mfpEnd,
transcriptionFactors=target.tf)
# tbl.fpAnnotated identifies 3 motifs in this one footprint, and each of them is associated with SATB2
# we therefore expect 3 rows, each describing the same footprint, and each with a different motif
checkEquals(dim(tbl), c(3, 18))
checkEquals(tbl$tfsMatched, rep("SATB2;CUX1", 3))
checkEquals(tbl$motifName, c("MA0679.1", "MA0756.1", "MA0757.1"))
} # test_1_footprint_two_tfs
#------------------------------------------------------------------------------------------------------------------------
test_emptyResults <- function()
{
printf("--- test_emptyResults")
# first test out a chromosomal region 1 base long
loc <- list(chr="chr1", start=3200, end=3200)
target.tf <- "ELF3" # mapped to the largest number of motifs in tbl.motifToMultipleGenes
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$start, endLoc=loc$end,
transcriptionFactors=target.tf)
checkEquals(nrow(tbl), 0)
# now repeat the SATB2 search (see above) with locs just 1 base too small, by incrementing start
loc <- tbl.fpAnnotated[grep("MA0679", tbl.fpAnnotated$motifName)[1], c("chr", "motifStart", "motifEnd")]
loc.orig <- loc
loc$motifStart <- loc$motifStart + 1
target.tf <- "SATB2"
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$motifStart, endLoc=loc$motifEnd,
transcriptionFactors=target.tf)
checkEquals(nrow(tbl), 0)
# now repeat the SATB2 search (see above) with locs just 1 base too small, by decrementing end
loc <- loc.orig
loc$motifEnd <- loc$motifEnd - 1
target.tf <- "SATB2"
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$start, endLoc=loc$end,
transcriptionFactors=target.tf)
checkEquals(nrow(tbl), 0)
} # test_emptyResults
#------------------------------------------------------------------------------------------------------------------------
test_360k_chr5_alzheimers_region <- function()
{
printf("--- test_360k_chr5_alzheimers_region")
loc <- list(chr="chr5", start=88640561, end=89001975)
target.tfs <- c("POU5F2", "SATB2", "HLF", "SOX4", "CUX2")
tbl <- findTFsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
chromosome=loc$chr, startLoc=loc$start, endLoc=loc$end,
transcriptionFactors=target.tfs)
checkEquals(dim(tbl), c(212, 18))
checkEquals(sort(unique(tbl$tfsMatched)), c("HLF", "POU5F2", "POU5F2;SOX4", "SATB2;CUX2", "SOX4"))
} # test_360k_chr5_alzheimers_region
#------------------------------------------------------------------------------------------------------------------------
test_findSNPsInFootprints <- function()
{
printf("--- test_findSNPsInFootprints")
snp.chromosome <- "chr5"
snp.loc.1 <- 88884587 # subset(tbl.snps, SNP %in% c("rs3814428", "rs770463", "rs446500"))$GRCh38.liftover
snp.loc.3 <- 88883758 # from 1000 genomes, should retrieve 6-motif footprint
tfs.all <- sort(unique(tbl.motifToMultipleGenes$tfs))
tfs <- c("ZNF263", "MAZ", "KLF16")
# search region
region.chromosome <- "chr5"
region.startLoc <- 88884570
region.endLoc <- 88884600 # chr5:88884570-88884600
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
region.chromosome, region.startLoc, region.endLoc,
snp.chromosome, snp.loc.1,
transcriptionFactors=tfs)
checkEquals(dim(tbl), c(7, 9))
checkEquals(sort(unique(tbl$tfsMatched)), c("KLF16", "MAZ", "ZNF263"))
# now retrieve a snp in a footprint, but notice that it misses the reported motif sequence
snp.loc.2 <- 88899133 # subset(tbl.snps, SNP %in% c("rs3814428", "rs770463", "rs446500"))$GRCh38.liftover
region.startLoc <- 88899100
region.endLoc <- 88899300
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
region.chromosome, region.startLoc, region.endLoc, # region of interest
snp.chromosome, snp.loc.2,
transcriptionFactors="TCF7")
checkEquals(dim(tbl), c(1,9))
x <- as.list(tbl[1,])
checkEquals(x$chr, "chr5")
checkEquals(x$mfpStart, 88899109)
checkEquals(x$mfpEnd, 88899140)
checkEquals(x$motifStart, 88899109)
checkEquals(x$motifEnd, 88899120)
checkEquals(x$sequence, "AAAGTTCAAAGC")
checkEquals(x$motifName, "MA0769.1")
checkEquals(x$snp, 88899133)
checkEquals(x$tfsMatched, "TCF7")
# now try combining the two snp.locs successfully tested above
tfs <- c("KLF16", "MAZ", "ZNF263", "TCF7")
region.startLoc <- 88884570
region.endLoc <- 88899300
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
region.chromosome, region.startLoc, region.endLoc, # region of interest
snp.chromosome, c(snp.loc.1, snp.loc.2),
transcriptionFactors=tfs)
checkEquals(dim(tbl), c(8, 9))
checkEquals(sort(unique(tbl$tfsMatched)), c("KLF16", "MAZ", "TCF7", "ZNF263"))
checkEquals(sort(unique(tbl$snp)), c(88884587, 88899133))
} # test_findSNPsInFootprints
#------------------------------------------------------------------------------------------------------------------------
# snp rs10038371 is just upstream chr5:88,903,407 of a footprint+motif at 88,903,378-403 (upstream from the perspective of the
# possibly regulated gene, MEF2C, which is encoded on the minus strand).
# here we look for that snp near that footprint, passing in a padding parameter which produces a liberal
# notion of overlapping
test_findSNPsNearFootprints <- function()
{
printf("--- test_findSNPsNearFootprints")
snp.chromosome <- "chr5"
snp.loc <- 88903407
target.region.chromosome <- "chr5"
target.region.startLoc <- 88903000
target.region.endLoc <- 88903500
tfs.all <- sort(unique(tbl.motifToMultipleGenes$tfs))
#tfs <- c("KLF16", "MAZ", "ZNF263", "TCF7")
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors=NA, # tfs.all,
padding=0)
checkEquals(nrow(tbl), 0)
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors=NA, # tfs.all,
padding=5)
checkEquals(dim(tbl), c(1,9))
# tfsMatched in the above table: ISL1;ISL2;LHX1;LHX2;LHX3;LHX4;LHX5;LHX6;LHX8;LHX9
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors="ISL1",
padding=5)
checkEquals(dim(tbl), c(1, 9))
# with very large padding, our snp rs10038371 "overlaps" with the footprint at 88,902,507:88,902534
# though of no likely biological interest, this helps to test the padding calculation
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc-1000, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors=NA, #tfs.all,
padding=1000)
checkEquals(dim(tbl), c(2, 9))
} # test_findSNPsNearFootprints
#------------------------------------------------------------------------------------------------------------------------
# a cryptic and transient bug, now no longer evident.
# originally, with all possible transcription factors considered, LHX1 was found (along with others in its gene family)
# for the region looked at here.
# but when LHX1 was considered alone, the results came back empty.
# this function ensures that they come back properly: one row describing an intersecting snp/fp on chr5 just upstream
# on the minus strand from mef2c
test_LXH1_matching <- function()
{
printf("--- test_LXH1_matching")
snp.chromosome <- "chr5"
snp.loc <- 88903407
target.region.chromosome <- "chr5"
target.region.startLoc <- 88903000
target.region.endLoc <- 88903500
tfs.all <- sort(unique(tbl.motifToMultipleGenes$tfs))
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors=NA, # tfs.all,
padding=5)
checkEquals(dim(tbl), c(1,9))
target.tf <- "LHX1"
checkTrue(grepl(target.tf, tbl[1, "tfsMatched"]))
tokens <- strsplit(tbl[1, "tfsMatched"], ";")[[1]]
checkTrue(match(target.tf, tokens) > 0)
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
target.region.chromosome, target.region.startLoc, target.region.endLoc,
snp.chromosome, snp.loc,
transcriptionFactors=target.tf,
padding=5)
checkEquals(dim(tbl), c(1,9))
} # test_LHX1_matching
#------------------------------------------------------------------------------------------------------------------------
explore_abca7_snps <- function()
{
printf("--- explore_abca7_snps")
f <- system.file(package="PrivateCoryData", "data", "abca7-associated-snps-1kg.tsv")
stopifnot(file.exists(f))
tbl.snps <- read.table(f, header=FALSE, sep="\t", as.is=TRUE)
colnames(tbl.snps) <- c("chrom", "loc")
displayBedTable(igv, tbl.snps[, c(1,2,2)], "1kg snps")
min.loc <- min(tbl.snps$loc)
max.loc <- max(tbl.snps$loc)
tbl.sub <- subset(tbl.fpAnnotated, chr=="chr19" & mfpStart >= (min.loc-500) & mfpEnd <= (max.loc + 1600))
displayBedTable(igv, tbl.sub[, c("chr", "mfpStart", "mfpEnd")], "footprints")
tbl <- findSNPsInFootprints(tbl.fpAnnotated, tbl.motifToMultipleGenes,
"chr19", min.loc, max.loc,
"chr19", sort(unique(tbl.snps$loc)),
transcriptionFactors=NA,
padding=10)
} # explore_abca7_snps
#------------------------------------------------------------------------------------------------------------------------
test_snpFootDisplay <- function()
{
printf("--- test_snpFootDisplay")
if(!exists("igv"))
igv <- igvR()
checkTrue(connected(igv))
#snpFootDisplay("chr5", start.loc, end.loc, genome="hg38")
vcfDirectory <- system.file(package="igvR", "extdata", "vcf")
vcfFile <- file.path(vcfDirectory, "chr22-sub.vcf.gz")
tbifile <- file.path(vcfDirectory, "chr22-sub.vcf.gz.tbi")
checkTrue(file.exists(vcfFile))
checkTrue(file.exists(tbiFile))
checkTrue(connected(igv))
samples.01 <- c("HG00096", "HG00097", "HG00099")
samples.02 <- c("HG00100", "HG00101")
start.loc <- 49902169
end.loc <- 49920831
goto(igv, "chr22", start.loc, end.loc)
displayVcfRegion(igv, "chr22", start.loc, end.loc, vcfFile, sampleIDs=samples.01)
displayVcfRegion(igv, "chr22", start.loc, end.loc, vcfFile, sampleIDs=samples.02)
displayVcfRegion(igv, "chr22", start.loc, end.loc, vcfFile)
} # test_snpFootDisplay
#------------------------------------------------------------------------------------------------------------------------
# "ITGA2B" is a problematic gene in /Volumes/local/Cory/for_Hongdong/trn.rtrim_isb_all.RData
test_tfGrabber <- function()
{
printf("--- test_tfGrabber")
dir <- system.file(package="snpFoot", "extdata")
checkTrue(file.exists(dir))
file <- "demoTRN.RData"
full.path <- file.path(dir, file)
checkTrue(file.exists(full.path))
var.names <- load(full.path, envir=.GlobalEnv)
var.name.trn.rtrim <- grep("^trn.rtrim", var.names, value=TRUE)
eval(parse(text=sprintf("trn.rtrim <<- %s", var.name.trn.rtrim[1])))
genes.of.interest <- intersect(rownames(trn.rtrim), c("MEF2C","ABCA7","CR1"))
promoterDistance <- 1000000
time.info <- system.time(x <- tfGrabber(genes.of.interest, trn.rtrim,
trnName="demo", promoterDist=promoterDistance))
print(time.info)
checkTrue(nrow(x$bed4igv) >= 10)
checkEquals(length(x$TRN), length(genes.of.interest))
} # test_tfGrabber
#------------------------------------------------------------------------------------------------------------------------
test_tfGrabber.new <- function()
{
printf("--- test_tfGrabber.new")
dir <- system.file(package="snpFoot", "extdata")
checkTrue(file.exists(dir))
file <- "demoTRN.RData"
full.path <- file.path(dir, file)
checkTrue(file.exists(full.path))
var.names <- load(full.path, envir=.GlobalEnv)
var.name.trn.rtrim <- grep("^trn.rtrim", var.names, value=TRUE)
eval(parse(text=sprintf("trn.rtrim <<- %s", var.name.trn.rtrim[1])))
genes.of.interest <- intersect(rownames(trn.rtrim), c("MEF2C","ABCA7","CR1"))
promoterDistance <- 1000000
x <- tfGrabber.new("MEF2C", trn.rtrim, trnName="demo", promoterDist=10)
checkEquals(dim(x), c(0,0))
x <- tfGrabber.new("MEF2C", trn.rtrim, trnName="demo", promoterDist=10000)
checkEquals(dim(x), c(4, 23))
x <- tfGrabber.new("MEF2C", trn.rtrim, trnName="demo", promoterDist=1000000)
checkEquals(dim(x), c(605, 23))
x <- tfGrabber.new("CR1", trn.rtrim, trnName="demo", promoterDist=1000000)
checkEquals(dim(x), c(0,0))
xx <- lapply(c("MEF2C", "ABCA7"), function(gene) tfGrabber.new(gene, trn.rtrim, trnName="demo", promoterDist=1e6))
xxx <- do.call("rbind", xx)
genes.44 <- c("ABCA7", "APOE", "BIN1", "CASS4", "CD2AP", "CELF1", "CLU",
"CR1", "EPHA1", "FERMT2", "HLA", "INPP5D", "MEF2C", "MS4A6A",
"NME8", "PICALM", "PTK2B", "SLC24A4/RIN3", "SORL1", "ZCWWPW1",
"ABI3", "IRF8", "MS4A14", "MS4A4A", "MS4A7", "PLCG1", "PLCG2",
"TGFB1", "TREM2", "TTC3", "UNC5C", "PILRB", "SLC24A4",
"FASLG", "TNFSF18", "PIEZO2", "NAPG", "APCDD1", "DNMBP",
"APP", "BACE1", "BACE2", "PSEN1", "PSEN2")
# 5 seconds elapsed
time.info <- system.time(fp.44.list <- lapply(genes.44, function(gene) tfGrabber.new(gene, trn.rtrim, trnName="demo", promoterDist=1e6)))
lapply(fp.44.list, length)
tbl.fp.44 <- do.call("rbind", fp.44.list)
dim(tbl.fp.44) # [1] 8470 23
table(tbl.fp.44$gene)
# ABCA7 ABI3 APCDD1 APOE BACE1 BACE2 BIN1 CASS4 CD2AP CELF1 DNMBP EPHA1
# 233 183 432 676 1117 57 18 62 167 177 57 30
# FERMT2 INPP5D MEF2C MS4A14 MS4A4A MS4A6A NAPG PIEZO2 PLCG1 PLCG2 PSEN1 PSEN2
# 411 185 605 122 16 49 445 181 121 158 453 177
# PTK2B SLC24A4 SORL1 TGFB1 TREM2 TTC3 UNC5C
# 300 84 218 1435 46 148 107
} # test_tfGrabber.new
#------------------------------------------------------------------------------------------------------------------------
test_run.tfGrabber.new <- function()
{
printf("--- test_run.tfGrabber.new")
load(system.file(package="PrivateCoryData", "extdata/trn10.genes44/genes.44.RData"))
checkTrue(exists("genes.44"))
checkEquals(length(genes.44), 44)
dir <- system.file(package="snpFoot", "extdata")
checkTrue(file.exists(dir))
file <- "demoTRN.RData"
full.path <- file.path(dir, file)
checkTrue(file.exists(full.path))
var.names <- load(full.path, envir=.GlobalEnv)
var.name.trn.rtrim <- grep("^trn.rtrim", var.names, value=TRUE)
eval(parse(text=sprintf("trn.rtrim <<- %s", var.name.trn.rtrim[1])))
checkEquals(dim(trn.rtrim), c(11822, 500))
genes.missing <- setdiff(genes.44, rownames(trn.rtrim))
genes.missing
# [1] "APP" "CLU" "CR1" "FASLG" "HLA" "IRF8" "MS4A7" "NME8" "PICALM" "PILRB" "TNFSF18" "ZCWWPW1"
intersect(genes.missing, colnames(trn.rtrim)) # [1] "IRF8"
genes.of.interest <- intersect(rownames(trn.rtrim), genes.44) # 32
sqlite.fileName <- file.path(getwd(), "grabber.sqlite")
if(file.exists(sqlite.fileName))
unlink(sqlite.fileName)
checkTrue(!file.exists(sqlite.fileName))
sql.tableName <- "tfFootprints"
tbl.fpDemo <- run.tfGrabber(genes.of.interest, trn.rtrim, "trnDemo", sqlite.fileName, sql.tableName)
db <- dbConnect(dbDriver("SQLite"), sqlite.fileName)
checkTrue(file.exists(sqlite.fileName))
result <- dbGetQuery(db, sprintf("select count(*) from %s", sql.tableName))
rowCount <- result[1,1]
checkTrue(rowCount > 8000)
checkTrue(rowCount < 9000)
tbl <- dbGetQuery(db, sprintf("select * from %s limit 3", sql.tableName))
checkEquals(dim(tbl), c(3, 24))
checkEquals(sort(colnames(tbl)), sort(c("TF", "beta", "chr", "end", "experimentCount", "fpEnd", "fpScore",
"fpStart", "gene", "id", "motifEnd", "motifFpOverlap", "motifName",
"motifScore", "motifStart", "motifStrand", "name", "promoterDist",
"pvalue", "qvalue", "row_names", "sequence", "start", "trnName")))
} # test_run.tfGrabber.new
#------------------------------------------------------------------------------------------------------------------------
test_intersectingLocs <- function()
{
printf("--- test_intersectingLocs")
chrom <- "chr1"
tbl.bed <- data.frame(chrom=chrom, start=1, end=10, stringsAsFactors=FALSE)
loc.good <- 6
loc.bad <- 16
checkEquals(intersectingLocs(tbl.bed, chrom, loc.good), loc.good)
checkEquals(intersectingLocs(tbl.bed, chrom, loc.bad), numeric(0))
checkEquals(intersectingLocs(tbl.bed, chrom, loc.bad, padding=10), loc.bad)
} # test_intersectingLocs
#------------------------------------------------------------------------------------------------------------------------
test_displayGWAS <- function()
{
printf("--- test_displayGWAS")
checkTrue(exists("tbl.gwas"))
if(!exists("igv"))
return(FALSE)
displayGWASTable(igv, tbl.gwas, "ADgwas2013")
} # test_displayGWAS
#------------------------------------------------------------------------------------------------------------------------
test_intersectMarietSnpsWithGWAS <- function()
{
printf("--- test_intersectMarietSnpsWithGWAS")
if(!exists("igv"))
return(FALSE)
# first, display all mariette snps
file <- system.file(package="PrivateCoryData", "extdata", "mef2c-related-snps-from-mariette.tsv")
checkTrue(file.exists(file))
tbl.mariette <- read.table(file, sep="\t", header=TRUE, as.is=TRUE)
mariette.snps <- tbl.mariette$BP
displaySnps(igv, "chr5", mariette.snps, "mariette.snps")
goto(igv, "chr5", 87850803, 88515362)
# now load all the footprints from 10trns, 44 genes
file <- system.file(package="PrivateCoryData", "extdata", "trn10.genes44", "44genes-10trns-bedTable.RData")
checkTrue(file.exists(file))
load(file) # it's called tbl.out, with almost 18k rows
displayBedTable(igv, tbl.out, "44genes.10trns") # igv will ask if you want to create an index. say yes
squishTrack(igv, "44genes.10trns.bed")
# create a track of intersection of gwas snps with 44genes.10trns footprints
gwas.chr5.snps <- subset(tbl.gwas, CHR=="chr5")$BP
gwas.chr5.snps.in.fp <- intersectingLocs(tbl.out, "chr5", gwas.chr5.snps, padding=10) # just two
displaySnps(igv, "chr5", gwas.chr5.snps.in.fp, "gwas.sps.in.fp")
# now do the same with mariette snps
mariette.snps.in.or.near.fp <- intersectingLocs(tbl.out, "chr5", mariette.snps, padding=50) # also just two
displaySnps(igv, "chr5", mariette.snps.in.or.near.fp, "mariette.sps.in.fp")
} # test_intersectMarietSnpsWithGWAS
#------------------------------------------------------------------------------------------------------------------------
test_displayAllEncodeFootprints <- function()
{
printf("--- test_displayAllEncodeFootprints")
if(!exists("igv"))
return()
checkTrue(exists("tbl.fpAnnotated"))
tbl.chr5 <- subset(tbl.fpAnnotated, chr=="chr5")
displayBedTable(igv, tbl.chr5[, c("chr", "mfpStart", "mfpEnd", "motifName")], "from.encode")
} # test_displayAllEncodeFootprints
#------------------------------------------------------------------------------------------------------------------------
test_runLift.hg19to38 <- function()
{
printf("--- test_runLift.hg19to38")
small.bed.file <- system.file(package="PrivateCoryData", "extdata/elizabethBlue/chr1_CU0039F.shared.bed")
expected.output.file <- "./chr1_CU0039F.shared.hg38.bed"
if(file.exists("expected.output.file"))
unlink(expected.output.file)
runLift.hg19to38(small.bed.file)
checkTrue(file.exists(expected.output.file))
tbl <- read.table(expected.output.file, sep="\t", as.is=TRUE)
checkEquals(dim(tbl), c(30, 4))
checkEquals(tbl[1,1], "chr1")
checkEquals(tbl[1,2], 172863016)
checkEquals(tbl[1,3], 172863016)
checkEquals(tbl[1,4], "1:172832156")
if(file.exists("expected.output.file"))
unlink(expected.output.file)
} # test_runLift.hg19to38
#------------------------------------------------------------------------------------------------------------------------
test_selectBestAmongDuplicates <- function()
{
printf("--- test_selectBestAmongDuplicates")
file.path <- system.file(package="snpFoot", "extdata", "footprints.tbl.withDups.RData")
load(file.path)
tbl.best <- selectBestAmongDuplicates(tbl.fpoi.snp)
checkEquals(nrow(tbl.best), 6)
} # test_selectBestAmongDuplicates
#------------------------------------------------------------------------------------------------------------------------
if(!interactive()) runTests()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.