Commit 77da31bb authored by Rauschenberger's avatar Rauschenberger
Browse files

automation

parent e78a60f1
......@@ -650,6 +650,9 @@ map.snps <- function(gene.chr,gene.start,gene.end,snp.chr,snp.pos,dist=10^3){
if(length(gene.chr)!=length(gene.start)|length(gene.chr)!=length(gene.end)){
stop("Invalid.",call.=FALSE)
}
if(!is.numeric(snp.chr)|!is.numeric(snp.pos)){
stop("Invalid.",call.=FALSE)
}
p <- length(gene.start)
x <- data.frame(from=integer(length=p),to=integer(length=p))
pb <- utils::txtProgressBar(min=0,max=p,style=3)
......@@ -794,9 +797,10 @@ test.single <- function(Y,X,map,i,limit=NULL,steps=NULL,rho=c(0,0.5,1)){
# extract data
ys <- map$exons[[i]]
y <- Y[,ys,drop=FALSE]; rm(Y)
y <- Y[,ys,drop=FALSE]
xs <- seq(from=map$snps$from[i],to=map$snps$to[i],by=1)
x <- X[,xs,drop=FALSE]; rm(X)
x <- X[,xs,drop=FALSE]
rm(Y,X); silent <- gc()
# test effects
pvalue <- rep(x=NA,times=length(rho))
......
This diff is collapsed.
......@@ -40,8 +40,9 @@ spliceQTL::get.exons.bbmri(path=path)
On the virtual machine, execute this chunk to test for alternative splicing. See the documentation of the R package spliceQTL for further information.
```{r analysis,eval=FALSE}
for(chr in 1:22){
for(chr in 22:1){
for(data in c("Geuvadis","LLS")){
rm(list=setdiff(ls(),c("data","chr","path"))); gc()
set.seed(1)
cat("Analysing",data,chr,":",as.character(Sys.time()),"\n")
......@@ -59,13 +60,13 @@ for(chr in 1:22){
exons <- list$exons; snps <- list$snps; rm(list)
cat("Adjusting samples:","\n")
#exons <- spliceQTL::adjust.samples(x=exons) # slow!
exons <- spliceQTL::adjust.samples(x=exons) # slow!
exons <- asinh(x=exons)
cat("Adjusting covariates:","\n")
names <- strsplit(x=colnames(exons),split="_") # exon names
length <- sapply(names,function(x) as.integer(x[[3]])-as.integer(x[[2]])) # exon length
#exons <- spliceQTL::adjust.covariates(x=exons,group=gene_id,offset=length) # slow!
exons <- spliceQTL::adjust.covariates(x=exons,group=gene_id,offset=length) # slow!
# subset chromosome
cond <- sapply(strsplit(x=colnames(exons),split="_"),function(x) x[[1]]==chr)
......@@ -79,8 +80,8 @@ for(chr in 1:22){
cat("Mapping SNPs:","\n")
names <- strsplit(x=colnames(snps),split=":")
snp.chr <- sapply(names,function(x) x[[1]])
snp.pos <- sapply(names,function(x) x[[2]])
snp.chr <- sapply(names,function(x) as.integer(x[[1]]))
snp.pos <- sapply(names,function(x) as.integer(x[[2]]))
map$snps <- spliceQTL::map.snps(gene.chr=map$genes$chr,
gene.start=map$genes$start,
gene.end=map$genes$end,
......@@ -90,8 +91,9 @@ for(chr in 1:22){
map <- spliceQTL::drop.trivial(map=map)
cat("Testing:",as.character(Sys.time())," -> ")
rm(list=setdiff(ls(),c("exons","snps","map","data","chr"))); gc()
pvalue <- spliceQTL::test.multiple(Y=exons,X=snps,map=map,rho=c(0,0.5),spec=1) # slow!
rm(list=setdiff(ls(),c("exons","snps","map","data","chr","path"))); gc()
# Use rho=c(0,0.5)
pvalue <- spliceQTL::test.multiple(Y=exons,X=snps,map=map,rho=0,spec=4) # slow!
cat(as.character(Sys.time()),"\n")
save(object=pvalue,file=file.path(path,paste0("pval.",data,".chr",chr,".RData")))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment