##Make a fasta file out of current probe.tab file #awk '{print ">" $2 "_" $3 "_" $4 "\n" $10}' MoGene-1_0-st-v1.probe.tab > MoGene_probes.fa ## create single fasta file from Mouse Genome #cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chrX.fa chrY.fa > ~/MGenome_filtered.fa ## create bowtie index for mm9 #nohup ~/software/bowtie-0.10.0/bowtie-build MGenome_filtered.fa ~/software/bowtie-0.10.0/indexes/mm9 & ##Run bowtie something like this: #bowtie --best --all -f -y -p 10 -v 0 -m 10 --unfa unaligned_mm9.fasta --maxfa multiple_mm9.fasta ~/software/bowtie-0.10.0/indexes/mm9 ./MoGene_probes.fa ./MoGene_mm9.map pr <- read.table("MoGene_mm9.map",sep="\t",header=FALSE,comment.char="",stringsAsFactors=FALSE) psid <- t(sapply(pr$V1, function(u) strsplit(u,"_",fixed=TRUE)[[1]], USE.NAMES=FALSE)) pr <- cbind(psid,pr[,-1],stringsAsFactors=FALSE) colnames(pr) <- c("transcript_cluster_id","probe.x","probe.y","strand","seqname","start") pr$probe.x <- as.numeric(pr$probe.x) pr$probe.y <- as.numeric(pr$probe.y) o <- order(pr$seqname, pr$start) pr <- pr[o,] inds <- split(seq_len(nrow(pr)), pr$transcript_cluster_id) keepInd <- lapply(inds, function(u) { tt <- table(pr$seqname[u],pr$strand[u]) w <- which(tt==max(tt), arr.ind=TRUE) # keep probes that map to the best chromosome keep <- pr$strand[u]==colnames(tt)[w[1,2]] & pr$seqname[u]==rownames(tt)[w[1,1]] u <- u[keep] # keep only probes that map once to the best chromosome ind <- pr$probe.x[u]*10000+pr$probe.y[u] tt1 <- table(ind) keep <- tt1[as.character(ind)]==1 u[keep] }) keepInd <- unlist(keepInd, use.names=FALSE) probetab <- pr[keepInd,] probetab$stop <- probetab$start+24 rm(pr,psid,inds,keepInd,o); gc() # save some space probetab$strand[probetab$strand=="+"] <- "P" probetab$strand[probetab$strand=="-"] <- "+" probetab$strand[probetab$strand=="P"] <- "-" write.table(probetab,file="MoGene_mm9_probetab.txt",row.names=FALSE,sep="\t",quote=FALSE) # create CDF from scratch source("http://bioinf.wehi.edu.au/folders/mrobinson/exon/flat2Cdf.R") flat2Cdf(file="MoGene_mm9_probetab.txt", chipType="MoGene-1_0-st-v1", rows=1050, cols=1050, xynames=c("probe.x","probe.y"), tags="mm9",col.class=c("integer","character")[c(2,1,1,2,2,1,1)], verbose=TRUE, sep="\t", ucol=1, gcol=1)