Mercurial > repos > greg > ideas_genome_tracks
changeset 0:a842ab7e3bba draft
Uploaded
| author | greg |
|---|---|
| date | Wed, 15 Nov 2017 13:00:07 -0500 |
| parents | |
| children | fe50d06f7052 |
| files | .shed.yml ideas_genome_tracks.R ideas_genome_tracks.xml test-data/e001-h3k4me3.bigwig test-data/test.bed |
| diffstat | 5 files changed, 298 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.shed.yml Wed Nov 15 13:00:07 2017 -0500 @@ -0,0 +1,11 @@ +name: ideas_genome_tracks +owner: greg +description: | + Contains a tool that creates genome track visualizations from IDEAS state files. +homepage_url: http://sites.stat.psu.edu/~yzz2/IDEAS/ +long_description: | + Contains a tool that creates genome track visualizations from IDEAS state files. +remote_repository_url: +type: unrestricted +categories: +- Epigenetics
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ideas_genome_tracks.R Wed Nov 15 13:00:07 2017 -0500 @@ -0,0 +1,174 @@ +#!/usr/bin/env Rscript + +suppressPackageStartupMessages(library("optparse")) + +option_list <- list( + make_option(c("-b", "--build"), action="store", dest="build", help="Genome build"), + make_option(c("-c", "--chrom_len_file"), action="store", dest="chrom_len_file", help="Chromosome length file"), + make_option(c("-d", "--header"), action="store", dest="header", default=NULL, help="Track header"), + make_option(c("-e", "--state_name"), action="store", dest="state_name", help="State name"), + make_option(c("-f", "--hub_id"), action="store", dest="hub_id", help="Not sure what this is"), + make_option(c("-i", "--cell_info"), action="store", dest="cell_info", default=NULL, help="Not sure what this is"), + make_option(c("-n", "--hub_name"), action="store", dest="hub_name", default=NULL, help="Not sure what this is"), + make_option(c("-p", "--input_dir_para"), action="store", dest="input_dir_para", help="Directory containing .para outputs from IDEAS"), + make_option(c("-q", "--input_dir_state"), action="store", dest="input_dir_state", help="Directory containing .state outputs from IDEAS"), + make_option(c("-u", "--output_track_db"), action="store", dest="output_track_db", help="Output track db file"), + make_option(c("-v", "--output_build"), action="store", dest="output_build", help="Output genome build file"), + make_option(c("-w", "--output_hub"), action="store", dest="output_hub", help="Output hub file") +) + +parser <- OptionParser(usage="%prog [options] file", option_list=option_list) +args <- parse_args(parser, positional_arguments=TRUE) +opt <- args$options + +create_color_scheme = function(input_dir_para) { + # Create the color scheme. + para_files <- list.files(path=input_dir_para, full.names=True); + mc = NULL; + x = read.table(opt$parameter_files, comment="!", nrows=1); + l = as.integer(regexpr("\\*", as.matrix(x))); + l = min(which(l>0)) - 2; + x = as.matrix(read.table(opt$parameter_files)); + if (length(para_files) > 1) { + for (i in 2:length(para_files)) { + x = x + as.matrix(read.table(paste(para_files[i], ".para", sep=""))); + } + } + p = x[,1] / sum(x[,1]); + m = array(as.matrix(x[,1:l+1] / x[,1]), dim=c(dim(x)[1], l)); + return get_state_color(m, mc); +} + +get_state_color = function(statemean, markcolor=NULL) +{ + if (length(markcolor) == 0) { + markcolor = rep("", dim(statemean)[2]); + markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:dim(statemean)[2]-1)/dim(statemean)[2], 1, 1); + markcolor = t(col2rgb(markcolor)); + } + rg = apply(statemean, 1, range); + mm = NULL; + for(i in 1:dim(statemean)[1]) { + mm = rbind(mm, (statemean[i,]-rg[1,i])/(rg[2,i]-rg[1,i]+1e-10)); + } + mm = mm^6; + mm = mm / (apply(mm, 1, sum) + 1e-10); + mycol = mm%*%markcolor; + s = apply(statemean, 1, max); + s = (s - min(s)) / (max(s) - min(s) + 1e-10); + h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])})); + h[,2] = h[,2] * s; + h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])}); + rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}), h); + return(rt); +} + +create_track = function(input_dir_state, chrom_len_file, outpref, state_color, header, state_name) { + state_files <- list.files(path=input_dir_state, full.names=True); + genome_size = read.table(chrom_len_file); + g = NULL; + for(i in state_files) { + message(paste(i, " ", sep=""), appendLF=F); + tg = as.matrix(fread(i)); + t = NULL; + for(j in 1:dim(genome_size)[1]) { + t = c(t, which(tg[,2]==as.character(genome_size[j, 1]) & as.numeric(tg[,4]) > as.numeric(genome_size[j, 2]))); + } + if (length(t) > 0) { + tg = tg[-t,]; + } + t = which(is.na(match(tg[,2], genome_size[,1]))==T); + if (length(t)>0) { + tg = tg[-t,]; + } + print(c(dim(g),dim(tg))); + g = rbind(g, tg); + } + message("Done"); + uchr = sort(unique(as.character(g[,2]))); + g1 = NULL; + for(i in uchr){ + t = which(g[,2]==i); + g1 = rbind(g1, g[t[order(as.integer(g[t, 3]))],]); + } + g = NULL; + chr = as.character(g1[,2]); + posst = as.numeric(g1[,3]); + posed = as.numeric(g1[,4]); + state = as.matrix(g1[,5:(dim(g1)[2]-1)]); + if (length(state_name)==0) { + state_name=0:max(state); + } + L = dim(g1)[1]; + n = dim(state)[2]; + if (length(header) > 0) { + colnames(g1) = header; + } + cells = as.character(colnames(g1)[5:(dim(g1)[2]-1)]); + g1 = NULL; + message("Generating tracks"); + options(scipen=999); + tt = which(chr[2:L]!=chr[2:L-1]); + tt = c(tt,which(posst[2:L]!=posed[2:L-1])); + tt = sort(unique(tt)); + for(i in 1:n) { + tstate = state[,i]; + t = c(tt,which(tstate[2:L]!=tstate[2:L-1])); + t = sort(unique(t)); + t0 = c(0, t) + 1; + t = c(t, L); + np = cbind(chr[t], posst[t0], posed[t], tstate[t]); + x = cbind(np[,1:3], state_name[as.integer(np[,4])+1], 1000, ".", np[,2:3], state_color[as.numeric(np[,4])+1]); + write.table(as.matrix(x), paste(outpref, i, "bed1", sep="."), quote=F, row.names=F, col.names=F); + print(x[1,]); + system(paste("bedToBigBed ", outpref, ".", i, ".bed1 ", chrom_len_file, " ", outpref, ".", i, ".bb", sep="")); + system(paste("rm ", paste(outpref, i, "bed1", sep="."))); + } + return(cells); +} + +create_track_db = function(input_dir_state, chrom_len_file, track_file_name, state_color, header, state_name, hub_name, cell_info) { + cells = create_track(input_dir_state, chrom_len_file, track_file_name, state_color, header, state_name); + if (length(cell_info) == 0) { + cell_info = cbind(cells, cells, cells, "#000000"); + cell_info = array(cell_info, dim=c(length(cells), 4)); + } + cell_info = as.matrix(cell_info); + track_db = NULL; + for(i in 1:length(cells)) { + ii = which(cells[i] == cell_info[,1]); + if (length(ii) == 0) { + next; + } + ii = ii[1]; + track_db = c(track_db, paste("track bigBed", i, sep="")); + track_db = c(track_db, paste("priority", ii)); + track_db = c(track_db, "type bigBed 9 ."); + track_db = c(track_db, "itemRgb on"); + track_db = c(track_db, "maxItems 100000"); + track_db = c(track_db, paste("shortLabel", cell_info[ii, 2])); + track_db = c(track_db, paste("longLabel", paste(hub_name, cell_info[ii, 3]))); + track_db = c(track_db, paste("color", paste(c(col2rgb(cell_info[ii, 4])), collapse=","))); + track_db = c(track_db, "visibility dense"); + track_db = c(track_db, ""); + } +} + +if (length(opt$hub_name) == 0) { + hub_name = opt$hub_id; +} + +# Create the tracks output directory. +tracks_output_dir = paste("tracks_", opt$hub_id, "/", sep=""); +dir.create(tracks_output_dir, showWarnings=FALSE); + +# Create the color scheme. +state_color <- create_color_scheme(opt$input_dir_para); + +# Create the tracks. +track_db <- create_track_db(opt$input_dir_state, opt$chrom_len_file, paste(tracks_output_dir, opt$hub_id, sep=""), state_color, opt$header, opt$state_name, opt$hub_name, opt$cell_info); + +# Write the outputs. +write.table(track_db, opt$output_track_db, quote=F, row.names=F, col.names=F); +write.table(c(paste("genome", opt$build), opt$output_build), paste(tracks_output_dir, "genomes_", opt$hub_id, ".txt", sep=""), quote=F, row.names=F, col.names=F); +write.table(c(paste("hub", opt$hub_id), paste("shortLabel", opt$hub_id), paste("longLabel", opt$hub_name), paste("genomesFile genomes_", opt$hub_id, ".txt", sep=""), "email yzz2 at psu.edu"), opt$output_hub, quote=F, row.names=F, col.names=F);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ideas_genome_tracks.xml Wed Nov 15 13:00:07 2017 -0500 @@ -0,0 +1,63 @@ +<tool id="ideas_genome_tracks" name="IDEAS genome tracks" version="1.0.0"> + <description></description> + <requirements> + <requirement type="package" version="2.4.27">bedops</requirement> + <requirement type="package" version="1.10.4">r-data.table</requirement> + <requirement type="package" version="1.4.4">r-optparse</requirement> + <requirement type="package" version="332">ucsc-bedtobigbed</requirement> + </requirements> + <command detect_errors="exit_code"><![CDATA[ +#import os +#set dbkey = None +#set input_dir_para = os.path.abspath('input_dir_para') +#set input_dir_state = os.path.abspath('input_dir_state') +mkdir $input_dir_para && +mkdir $input_dir_state && +#for $i in $input: + #if $dbkey is None: + #set dbkey = $i.metadata.dbkey + #end if + #set filename = $i.file_name + #set name = $i.name + #if $name.endswith(".para"): + ln -s $filename $input_dir_para/$name && + #else if $name.endswith(".state"): + n -s $filename $input_dir_state/$name && + #end if +#end for +Rscript '$__tool_directory__/ideas_genome_tracks.R' +-b $dbkey +-c '$chromInfo' +-d '$track_header' +-f test +-p '$input_dir_para' +-q '$input_dir_state' +-u '$output_track_db' +-v '$output_build' +-w '$output_hub' + ]]></command> + <inputs> + <param name="input" format="fasta" type="data_collection" collection_type="list" label="IDEAS files" /> + <param name="track_header" type="text" value="" label="Track header"/> + </inputs> + <outputs> + <data name="output_track_db" format="bigbed" label="${tool.name} (track db) on ${on_string}"/> + <data name="output_build" format="txt" label="${tool.name} (genome build) on ${on_string}"/> + <data name="output_hub" format="txt" label="${tool.name} (hub) on ${on_string}"/> + </outputs> + <tests> + <test> + </test> + </tests> + <help> +**What it does** + +----- + +**Required options** + + </help> + <citations> + <citation type="doi">10.1093/nar/gkw278</citation> + </citations> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/test.bed Wed Nov 15 13:00:07 2017 -0500 @@ -0,0 +1,50 @@ +chr1 21819600 21819800 R100001 +chr1 21819800 21820000 R100002 +chr1 21820000 21820200 R100003 +chr1 21820200 21820400 R100004 +chr1 21820400 21820600 R100005 +chr1 21820600 21820800 R100006 +chr1 21820800 21821000 R100007 +chr1 21821000 21821200 R100008 +chr1 21821200 21821400 R100009 +chr1 21821400 21821600 R100010 +chr1 21821600 21821800 R100011 +chr1 21821800 21822000 R100012 +chr1 21822000 21822200 R100013 +chr1 21822200 21822400 R100014 +chr1 21822400 21822600 R100015 +chr1 21822600 21822800 R100016 +chr1 21822800 21823000 R100017 +chr1 21823000 21823200 R100018 +chr1 21823200 21823400 R100019 +chr1 21823400 21823600 R100020 +chr1 21823600 21823800 R100021 +chr1 21823800 21824000 R100022 +chr1 21824000 21824200 R100023 +chr1 21824200 21824400 R100024 +chr1 21824400 21824600 R100025 +chr1 21824600 21824800 R100026 +chr1 21824800 21825000 R100027 +chr1 21825000 21825200 R100028 +chr1 21825200 21825400 R100029 +chr1 21825400 21825600 R100030 +chr1 21825600 21825800 R100031 +chr1 21825800 21826000 R100032 +chr1 21826000 21826200 R100033 +chr1 21826200 21826400 R100034 +chr1 21826400 21826600 R100035 +chr1 21826600 21826800 R100036 +chr1 21826800 21827000 R100037 +chr1 21827000 21827200 R100038 +chr1 21827200 21827400 R100039 +chr1 21827400 21827600 R100040 +chr1 21827600 21827800 R100041 +chr1 21827800 21828000 R100042 +chr1 21828000 21828200 R100043 +chr1 21828200 21828400 R100044 +chr1 21828400 21828600 R100045 +chr1 21829000 21829200 R100046 +chr1 21829400 21829600 R100047 +chr1 21829600 21829800 R100048 +chr1 21829800 21830000 R100049 +chr1 21830000 21830200 R100050