26 September 2011

PostScript as a Programming Language for Bioinformatics: mynotebook

"PostScript (PS) is an interpreted, stack-based programming language. It is best known for its use as a page description language in the electronic and desktop publishing areas."[wikipedia]. In this post, I'll show how I've used to create a simple and lightweight view of the genome.

Introduction: just a simple postscript program

The following PS program fills a rectangular gray shape; You can display the result using ghostview, a2ps, etc... 
%!PS
newpath
50 50 moveto
0 100 rlineto
100 0 rlineto
0 -100 rlineto
closepath
0.5 setgray
fill
showpage

Some global variables

The page width

/screenWidth 1000 def

The page width

/screenHeight 1000 def

The minimum 5' position

/minChromStart 1E9 def

The maximum 3' position

/maxChromEnd -1 def

The size of a genomic feature

/featureHeight 20 def

The distance between two 'ticks' for drawing the orientation

/ticksx 20 def

The font size

/theFontSize 9 def
The variable knownGene is a PS array of genes.

/knownGene [
[(uc002zkr.3) (chr22) (-) 161242...
...]
] def

Each Gene is a PS array holding the structure of the UCSC knownGene table, that is to say: name , chromosome, txStart, txEnd, cdsStart, cdsEnd, exonStarts, exonEnds:

[(uc002zmh.2) (chr22) (-) 17618410 17646177 17618910 17646134
   [17618410 17619439 17621948 17623987 17625913 17629337 17630431 17646098 ]
   [17619247 17619628 17622123 17624021 17626007 17629450 17630635 17646177 ]
]
. a simple command line can be used to fetch those data:
%  curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/knownGene.txt.gz" |\
gunzip -c | grep chr22 | head -n 20 |\
awk '{printf("[(%s) (%s) (%s) %s %s %s %s [%s] [%s] ]\n",$1,$2,$3,$4,$5,$6,$7,$9,$10);}' |\
tr "," " " > result.txt

Some utilities

converting a PS object to string

/toString
{
20 string cvs 
} bind def

Converting a string to interger (loop over each character an increase the current value)

/toInteger
{
3 dict begin
/s exch def
/i 0 def
/n 0 def
s {
  n 10 mul /n exch def
  s i get 48 sub n add /n exch def %48=ascii('0')
  i 1 add /i exch def
  } forall
n % leave n on the stack
end
} bind def

Convert a genomic position to a index on the page 'x' axis

/convertPos2pixel
{
minChromStart sub maxChromEnd minChromStart sub div screenWidth mul
} bind def

Extract the chromosome (that is to say, extract the 1st element of the current array on the stack)

/getChrom
{
1 get
} bind def

Create a hyperlink to the UCSC genome browser

/getHyperLink
{
3 dict begin
/E exch def %% END
/S exch def %% START
/C exch def %% CHROMOSOME
[ (http://genome.ucsc.edu/cgi-bin/hgTracks?position=) C (:) S toString (-) E toString (&) (&db=hg19) ] concatstringarray
end
} bind def

Paint a rectangle

/box
{
4 dict begin
/height exch def
/width exch def
/y exch def
/x exch def
x y moveto
width 0 rlineto
0 height rlineto
width -1 mul 0 rlineto
0 height -1 mul rlineto
end
} bind def

Paint a gray gradient

/gradient
{
4 dict begin
/height exch def
/width exch def
/y exch def
/x exch def
/i 0 def
height 2 div /i exch def

0 1 height 2 div {
	1 i height 2.0 div div sub setgray
	newpath
	x  
	y height 2 div i sub  add
	width
	i 2 mul
	box
	closepath
	fill
	i 1 sub /i exch def
	}for
newpath
0 setgray
0.4 setlinewidth
x y width height box
closepath
stroke
end
} bind def

Methods extracting a data about the current gene on the PS stack.

Extract the transcription start:

/getTxStart
{
3 get
} bind def

Extract the transcription end:

/getTxEnd
{
4 get
} bind def

Extract the CDS start:

/getCdsStart
{
5 get
} bind def

Extract the transcription end:

/getCdsEnd
{
6 get
} bind def

Extract the strand:

/getStrand
{
2 get (+) eq {1} {-1} ifelse
} bind def
Get the gene name

/getKgName
{
0 get
} bind def

Get the number of exons:

/getExonCount
{
7 get length
} bind def

Get the start position of the i-th exon:

/getExonStart
{
2 dict begin
/i exch def
/gene exch def
gene 7 get i get
end
} bind def

Get the end position of the i-th exon:

/getExonEnd
{
2 dict begin
/i exch def
/gene exch def
gene 8 get i get
end
} bind def

Should we draw this gene on the page ?

/isVisible
{
1 dict begin
/gene exch def
minChromStart gene getTxEnd gt 
	{
	false
	}
	{
	gene getTxStart maxChromEnd gt
		{
		false
		}
		{
		true
		}ifelse
	}ifelse
end
}bind def

Methods for an array of genes

Loop over the genes and extract the lowest 5' index:

/getMinChromStart
{
3 dict begin
/genes exch def
/pos 10E9 def
/i 0 def
genes length {
	genes i get getTxStart pos min /pos  exch def
	i 1 add /i exch def
	}repeat
pos
end
} bind def

Loop over the genes and extract the highest 3' index:

/getMaxChromEnd
{
3 dict begin
/genes exch def
/pos -1E9 def
/i 0 def
genes length {
	genes i get getTxEnd pos max /pos  exch def
	i 1 add /i exch def
	}repeat
pos
end
} bind def

Painting ONE Gene

/paintGene
{
5 dict begin
/gene exch def %% the GENE argument
/midy featureHeight 2.0 div def %the middle of the row
/x0 gene getTxStart convertPos2pixel def % 5' side of the gene in pixel
/x1 gene getTxEnd convertPos2pixel def % 3' side of the gene in pixel
/i 0 def
0.1 setlinewidth

1 0 0 setrgbcolor

newpath
x0 midy moveto
x1 midy lineto
closepath
stroke


% paint ticks
0 1 x1 x0 sub ticksx div{
	newpath
	gene getStrand 1 eq 
		{
		x0 ticksHeight sub i add midy ticksHeight add moveto
		x0 i add midy lineto
		x0 ticksHeight sub i add midy ticksHeight sub lineto
		}
	%else
		{
		x0 ticksHeight add i add midy ticksHeight add moveto
		x0 i add midy lineto
		x0 ticksHeight add i add midy ticksHeight sub lineto
		} ifelse
	stroke
	i ticksx add /i exch def
	} for

%paint Transcript start-end
0 0 1 setrgbcolor
newpath
gene getCdsStart convertPos2pixel
midy cdsHeight 2 div sub
gene getCdsEnd convertPos2pixel gene getCdsStart convertPos2pixel sub 
cdsHeight box
closepath
fill

% loop over exons
0 /i exch def
gene getExonCount
	{
	gene i getExonStart convertPos2pixel
	midy exonHeight 2 div sub
	gene i getExonEnd convertPos2pixel gene i getExonStart convertPos2pixel sub
	exonHeight gradient
	i 1 add /i exch def
	} repeat
0 setgray
gene getTxEnd convertPos2pixel 10 add midy moveto
gene getKgName show

%URL 
[ /Rect [x0 0 x1 1 add featureHeight]
/Border [1 0 0]
/Color [1 0 0]
/Action << /Subtype /URI /URI gene getChrom gene getTxStart gene getTxEnd getHyperLink  >>
/Subtype /Link
/ANN pdfmark

end
} bind def

Paint all Genes

/paintGenes
{
3 dict begin
/genes exch def %the GENE argument (an array)
/i 0 def % loop iterator
/j 0 def % row iterator


% draw 10 vertical lines
i 0 /i exch def
0 setgray
0 1 10 {
	%draw a vertical line
	screenWidth 10 div i mul 0 moveto
	screenWidth 10 div i mul screenHeight lineto
	stroke
	% print the position at the top rotate by 90°
	screenWidth 10 div i mul 10 add screenHeight 5 sub moveto
	-90 rotate
	maxChromEnd minChromStart sub i 10 div mul minChromStart add toString show
	90 rotate
	i 1 add /i exch def
	} for

0 /i exch def
genes length {
	genes i get isVisible
		{
		gsave
		0 j  featureHeight 2 add mul translate
		genes i get paintGene
		j 1 add /j exch def
		grestore
		} if
	i 1 add /i exch def
	}repeat
end
} bind def

All in one: the postscript code

%!PS
%%Creator: Pierre Lindenbaum PhD plindenbaum@yahoo.fr http://plindenbaum.blogspot.com
%%Title: Genome Browser
%%CreationDate: 2011-09-25
%%BoundingBox: 0 0 1000 1000
%%Pages! 1
%
% curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/knownGene.txt.gz" | gunzip -c | head -n 20 | awk '{printf("[(%s) (%s) (%s) %s %s %s %s [%s] [%s] ]\n",$1,$2,$3,$4,$5,$6,$7,$9,$10);}' | tr "," " "
%
/pdfmark where {pop} {userdict /pdfmark /cleartomark load put} ifelse
/screenWidth 1000 def
/screenHeight 1000 def
/minChromStart 1E9 def
/maxChromEnd -1 def
/featureHeight 20 def
/ticksx 20 def
/theFontSize 9 def
/Courier findfont theFontSize scalefont setfont
%% http://en.wikibooks.org/wiki/PostScript_FAQ#How_to_concatenate_strings.3F
/concatstringarray % [(a) (b) ... (z)] --> (ab...z)
{ 0 1 index { length add } forall string
0 3 2 roll
{ 3 copy putinterval
length add
}
forall pop
} bind def
/toString
{
20 string cvs
} bind def
/toInteger
{
3 dict begin
/s exch def
/i 0 def
/n 0 def
s {
n 10 mul /n exch def
s i get 48 sub n add /n exch def %48=ascii('0')
i 1 add /i exch def
} forall
n
end
} bind def
/convertPos2pixel
{
minChromStart sub maxChromEnd minChromStart sub div screenWidth mul
} bind def
/getChrom
{
1 get
} bind def
/getTxStart
{
3 get
} bind def
/getTxEnd
{
4 get
} bind def
/getCdsStart
{
5 get
} bind def
/getCdsEnd
{
6 get
} bind def
/getStrand
{
2 get (+) eq {1} {-1} ifelse
} bind def
/getKgName
{
0 get
} bind def
/getExonCount
{
7 get length
} bind def
/getExonStart
{
2 dict begin
/i exch def
/gene exch def
gene 7 get i get
end
} bind def
/getExonEnd
{
2 dict begin
/i exch def
/gene exch def
gene 8 get i get
end
} bind def
/isVisible
{
1 dict begin
/gene exch def
minChromStart gene getTxEnd gt
{
false
}
{
gene getTxStart maxChromEnd gt
{
false
}
{
true
}ifelse
}ifelse
end
}bind def
/getHyperLink
{
3 dict begin
/E exch def
/S exch def
/C exch def
[ (http://genome.ucsc.edu/cgi-bin/hgTracks?position=) C (:) S toString (-) E toString (&) (&db=hg19) ] concatstringarray
end
} bind def
/getMinChromStart
{
3 dict begin
/genes exch def
/pos 10E9 def
/i 0 def
genes length {
genes i get getTxStart pos min /pos exch def
i 1 add /i exch def
}repeat
pos
end
} bind def
/getMaxChromEnd
{
3 dict begin
/genes exch def
/pos -1E9 def
/i 0 def
genes length {
genes i get getTxEnd pos max /pos exch def
i 1 add /i exch def
}repeat
pos
end
} bind def
/box
{
4 dict begin
/height exch def
/width exch def
/y exch def
/x exch def
x y moveto
width 0 rlineto
0 height rlineto
width -1 mul 0 rlineto
0 height -1 mul rlineto
end
} bind def
/gradient
{
4 dict begin
/height exch def
/width exch def
/y exch def
/x exch def
/i 0 def
height 2 div /i exch def
0 1 height 2 div {
1 i height 2.0 div div sub setgray
newpath
x
y height 2 div i sub add
width
i 2 mul
box
closepath
fill
i 1 sub /i exch def
}for
newpath
0 setgray
0.4 setlinewidth
x y width height box
closepath
stroke
end
} bind def
/ticksHeight
{
featureHeight 0.2 mul
} bind def
/cdsHeight
{
featureHeight 0.5 mul
} bind def
/exonHeight
{
featureHeight 0.8 mul
} bind def
%*********************************
%
% paint one gene
% @arg one knownGene array
% @return void
%
/paintGene
{
5 dict begin
/gene exch def %% the GENE argument
/midy featureHeight 2.0 div def %the middle of the row
/x0 gene getTxStart convertPos2pixel def % 5' side of the gene in pixel
/x1 gene getTxEnd convertPos2pixel def % 3' side of the gene in pixel
/i 0 def
0.1 setlinewidth
1 0 0 setrgbcolor
newpath
x0 midy moveto
x1 midy lineto
closepath
stroke
% paint ticks
0 1 x1 x0 sub ticksx div{
newpath
gene getStrand 1 eq
{
x0 ticksHeight sub i add midy ticksHeight add moveto
x0 i add midy lineto
x0 ticksHeight sub i add midy ticksHeight sub lineto
}
%else
{
x0 ticksHeight add i add midy ticksHeight add moveto
x0 i add midy lineto
x0 ticksHeight add i add midy ticksHeight sub lineto
} ifelse
stroke
i ticksx add /i exch def
} for
%paint Transcript start-end
0 0 1 setrgbcolor
newpath
gene getCdsStart convertPos2pixel
midy cdsHeight 2 div sub
gene getCdsEnd convertPos2pixel gene getCdsStart convertPos2pixel sub
cdsHeight box
closepath
fill
% loop over exons
0 /i exch def
gene getExonCount
{
gene i getExonStart convertPos2pixel
midy exonHeight 2 div sub
gene i getExonEnd convertPos2pixel gene i getExonStart convertPos2pixel sub
exonHeight gradient
i 1 add /i exch def
} repeat
0 setgray
gene getTxEnd convertPos2pixel 10 add midy moveto
gene getKgName show
%URL
[ /Rect [x0 0 x1 1 add featureHeight]
/Border [1 0 0]
/Color [1 0 0]
/Action << /Subtype /URI /URI gene getChrom gene getTxStart gene getTxEnd getHyperLink >>
/Subtype /Link
/ANN pdfmark
end
} bind def
%********************************************
%
% count the number of visible genes
%
%
/countVisibleGenes
{
3 dict begin
/genes exch def %the GENE argument (an array)
/i 0 def % loop iterator
/n 0 def % the count
genes length {
genes i get isVisible
{
n 1 add /n exch def
} if
i 1 add /i exch def
}repeat
n
end
} bind def
%********************************************
%
% draw an array of genes
%
% @arg an array of knownGene
% @return void
%
/paintGenes
{
3 dict begin
/genes exch def %the GENE argument (an array)
/i 0 def % loop iterator
/j 0 def % row iterator
% draw 10 vertical lines
i 0 /i exch def
0 setgray
0 1 10 {
%draw a vertical line
screenWidth 10 div i mul 0 moveto
screenWidth 10 div i mul screenHeight lineto
stroke
% print the position at the top rotate by 90°
screenWidth 10 div i mul 10 add screenHeight 5 sub moveto
-90 rotate
maxChromEnd minChromStart sub i 10 div mul minChromStart add toString show
90 rotate
i 1 add /i exch def
} for
0 /i exch def
genes length {
genes i get isVisible
{
gsave
0 j featureHeight 2 add mul translate
genes i get paintGene
j 1 add /j exch def
grestore
} if
i 1 add /i exch def
}repeat
end
} bind def
/knownGene [
[(uc002zkr.3) (chr22) (-) 16124263 16193004 16124263 16124263 [16124263 16162396 16186810 16187164 16192905 ] [16124973 16162487 16186946 16187302 16193004 ] ]
[(uc002zks.3) (chr22) (-) 16150261 16193004 16150261 16150261 [16150261 16162396 16186810 16187164 16189031 16189263 16190680 16192905 ] [16151821 16162487 16186946 16187302 16189143 16189411 16190791 16193004 ] ]
[(uc002zkt.2) (chr22) (+) 16162065 16172264 16162065 16162065 [16162065 16164481 16171951 ] [16162388 16164569 16172264 ] ]
[(uc002zku.2) (chr22) (-) 16179619 16181004 16179619 16179619 [16179619 ] [16181004 ] ]
[(uc002zkv.3) (chr22) (-) 16187164 16193004 16187164 16187164 [16187164 16189031 16189263 16190680 16192905 ] [16187302 16189143 16189378 16190791 16193004 ] ]
[(uc002zkw.2) (chr22) (-) 16240244 16240281 16240244 16240244 [16240244 ] [16240281 ] ]
[(uc002zkx.1) (chr22) (-) 16240300 16240339 16240300 16240300 [16240300 ] [16240339 ] ]
[(uc002zky.1) (chr22) (+) 16241086 16241125 16241086 16241086 [16241086 ] [16241125 ] ]
[(uc002zkz.1) (chr22) (-) 16242000 16242030 16242000 16242000 [16242000 ] [16242030 ] ]
[(uc002zla.1) (chr22) (-) 16243380 16243414 16243380 16243380 [16243380 ] [16243414 ] ]
[(uc002zlb.2) (chr22) (-) 16243908 16243948 16243908 16243908 [16243908 ] [16243948 ] ]
[(uc002zlc.1) (chr22) (-) 16245151 16245185 16245151 16245151 [16245151 ] [16245185 ] ]
[(uc002zld.2) (chr22) (-) 16245679 16245719 16245679 16245679 [16245679 ] [16245719 ] ]
[(uc002zle.1) (chr22) (-) 16248998 16249023 16248998 16248998 [16248998 ] [16249023 ] ]
[(uc002zlf.1) (chr22) (-) 16251234 16254941 16251234 16251234 [16251234 ] [16254941 ] ]
[(uc002zlg.1) (chr22) (-) 16256331 16287425 16256331 16256331 [16256331 16258184 16262903 16266928 16268136 16269872 16275206 16277747 16279194 16282144 16282477 16287253 ] [16256677 16258303 16262952 16267095 16268181 16269943 16275277 16277885 16279301 16282318 16282592 16287425 ] ]
[(uc002zlh.1) (chr22) (-) 16256331 16287425 16258185 16280411 [16256331 16258184 16266928 16268136 16269872 16275206 16277747 16279194 16280333 16282144 16282477 16287253 ] [16256677 16258303 16267095 16268181 16269943 16275277 16277885 16279301 16280589 16282318 16282592 16287425 ] ]
[(uc010gqp.2) (chr22) (-) 16256331 16287937 16258185 16287885 [16256331 16258184 16266928 16268136 16269872 16275206 16277747 16279194 16282144 16282477 16287253 ] [16256677 16258303 16267095 16268181 16269943 16275277 16277885 16279301 16282318 16282592 16287937 ] ]
[(uc002zlj.1) (chr22) (-) 16266928 16287937 16266930 16287390 [16266928 16268136 16269872 16275206 16277747 16279194 16282144 16282477 16287253 16287537 ] [16267095 16268181 16269943 16275277 16277885 16279301 16282318 16282592 16287425 16287937 ] ]
[(uc002zlk.2) (chr22) (+) 16274557 16278598 16274557 16274557 [16274557 16276480 ] [16275003 16278598 ] ]
[(uc002zll.1) (chr22) (+) 16373080 16377057 16373080 16373080 [16373080 16373829 16375448 ] [16373121 16373911 16377057 ] ]
[(uc011agd.1) (chr22) (-) 16448825 16449804 16448826 16449804 [16448825 ] [16449804 ] ]
[(uc002zln.1) (chr22) (+) 16492811 16492932 16492811 16492811 [16492811 ] [16492932 ] ]
[(uc011age.1) (chr22) (+) 17029052 17029078 17029052 17029052 [17029052 ] [17029078 ] ]
[(uc002zlo.1) (chr22) (+) 17029615 17029643 17029615 17029615 [17029615 ] [17029643 ] ]
[(uc002zlp.1) (chr22) (-) 17071647 17073700 17071766 17073440 [17071647 ] [17073700 ] ]
[(uc010gqq.2) (chr22) (+) 17082800 17095998 17082800 17082800 [17082800 17092547 17094966 17095588 ] [17083105 17092783 17095068 17095998 ] ]
[(uc002zlq.3) (chr22) (+) 17082800 17095998 17082800 17082800 [17082800 17092547 17094966 ] [17083105 17092783 17095998 ] ]
[(uc002zlr.2) (chr22) (+) 17082800 17129719 17082800 17082800 [17082800 17092547 17094966 17103730 17117929 17119468 17128057 17128552 17129416 ] [17083105 17092783 17095068 17103787 17117980 17119630 17128147 17128675 17129719 ] ]
[(uc002zls.1) (chr22) (+) 17082800 17179521 17082800 17082800 [17082800 17119468 17178385 ] [17083105 17119630 17179521 ] ]
[(uc002zlt.2) (chr22) (+) 17100506 17134580 17100506 17100506 [17100506 17100941 17103730 17117929 17119468 17128494 17129416 17134399 ] [17100610 17101050 17103787 17117980 17119630 17128675 17129538 17134580 ] ]
[(uc002zlu.2) (chr22) (-) 17227759 17229328 17227759 17227759 [17227759 17229165 ] [17228629 17229328 ] ]
[(uc002zlv.2) (chr22) (-) 17264312 17302584 17264508 17288963 [17264312 17280660 17288628 17302496 ] [17265299 17280914 17288973 17302584 ] ]
[(uc011agf.1) (chr22) (-) 17264312 17302584 17264508 17288963 [17264312 17280660 17288628 17302496 ] [17265299 17280914 17288976 17302584 ] ]
[(uc010gqr.1) (chr22) (+) 17308363 17310225 17308363 17308363 [17308363 17309431 ] [17308950 17310225 ] ]
[(uc011agg.1) (chr22) (-) 17385314 17385395 17385314 17385314 [17385314 ] [17385395 ] ]
[(uc002zlw.2) (chr22) (-) 17442828 17489112 17443622 17489004 [17442828 17444614 17445655 17446067 17446989 17449187 17450832 17468849 17472762 17488830 ] [17443766 17444719 17445752 17446158 17447254 17449273 17451083 17469057 17473066 17489112 ] ]
[(uc010gqs.1) (chr22) (-) 17446991 17489112 17449218 17489004 [17446991 17449187 17468849 17472762 17488830 ] [17448371 17449273 17469006 17473066 17489112 ] ]
[(uc002zlx.1) (chr22) (+) 17517459 17539682 17517459 17517459 [17517459 17525762 17528213 17537962 ] [17518234 17525890 17528316 17539682 ] ]
[(uc002zly.2) (chr22) (+) 17565848 17591387 17565981 17590710 [17565848 17577951 17578686 17579664 17581244 17582885 17583028 17584383 17585615 17586480 17586742 17588616 17589196 ] [17566119 17577976 17578833 17579777 17581371 17582933 17583192 17584467 17585700 17586492 17586844 17588658 17591387 ] ]
[(uc010gqt.2) (chr22) (+) 17565848 17591387 17565981 17590710 [17565848 17577951 17578686 17579664 17581244 17582885 17583028 17584383 17585615 17589196 ] [17566119 17577976 17578833 17579777 17581371 17582933 17583192 17584467 17585700 17591387 ] ]
[(uc002zlz.1) (chr22) (+) 17593917 17596582 17593917 17593917 [17593917 ] [17596582 ] ]
[(uc002zmb.2) (chr22) (-) 17597188 17602213 17600280 17602017 [17597188 ] [17602213 ] ]
[(uc002zma.2) (chr22) (-) 17597188 17602257 17600280 17600952 [17597188 17602141 ] [17601033 17602257 ] ]
[(uc002zmc.2) (chr22) (+) 17602484 17612993 17602484 17602484 [17602484 17605544 17611251 17612504 ] [17602929 17605661 17611374 17612993 ] ]
[(uc002zmd.2) (chr22) (-) 17618410 17622467 17618910 17622286 [17618410 17619439 17621948 17622282 ] [17619247 17619628 17622123 17622467 ] ]
[(uc002zme.2) (chr22) (-) 17618410 17629450 17618910 17622070 [17618410 17619439 17621948 17625913 17629337 ] [17619247 17619628 17622123 17626007 17629450 ] ]
[(uc002zmf.2) (chr22) (-) 17618410 17640169 17618910 17640141 [17618410 17619439 17621948 17623987 17625913 17629337 17630431 17640015 ] [17619247 17619628 17622123 17624021 17626007 17629450 17630635 17640169 ] ]
[(uc002zmg.2) (chr22) (-) 17618410 17640169 17618910 17640141 [17618410 17621948 17623987 17640015 ] [17619247 17622123 17624021 17640169 ] ]
[(uc002zmh.2) (chr22) (-) 17618410 17646177 17618910 17646134 [17618410 17619439 17621948 17623987 17625913 17629337 17630431 17646098 ] [17619247 17619628 17622123 17624021 17626007 17629450 17630635 17646177 ] ]
] def
knownGene getMinChromStart 10 sub /minChromStart exch def
knownGene getMaxChromEnd 10 add /maxChromEnd exch def
systemdict /userChromStart known {
userChromStart toInteger /minChromStart exch def
} if
systemdict /userChromEnd known {
userChromEnd toInteger /maxChromEnd exch def
} if
screenHeight knownGene countVisibleGenes 1 add div 20.0 min /featureHeight exch def
0 0 screenWidth 1 sub screenHeight 1 sub box
0.95 setgray
fill
knownGene paintGenes
0 0 screenWidth screenHeight box
0 setgray
stroke
showpage
view raw genomebrowser.ps.txt hosted with ❤ by GitHub

Open the PS file in ghostview, evince, ...

Zooming ? Yes we can.

Ghostview has an option -Sname=string
       -Sname=string
       -sname=string
              Define  a  name  in  "systemdict"  with a given string as value.
              This is different from -d.

In my postscript file, the default values for minChromStart and maxChromEnd are overridden by the user's parameters:

systemdict /userChromStart known {
	userChromStart toInteger /minChromStart  exch def
	} if

systemdict /userChromEnd known {
	userChromEnd toInteger /maxChromEnd  exch def
	} if
That's it,

Pierre

Publié par Pierre Lindenbaum at 11:02 PM

Libellés : , , , , ,

3 comments:

Tim said...

That's really lightweight, great! Is there a way to pass the gene list in the command-line instead of copy-paste it in the ".ps" file?

Tuesday, 27 September, 2011
Pierre Lindenbaum said...

@Tim I don't know, I'm not a PS guru. I briefly searched on the web but found no answer.

Tuesday, 27 September, 2011
Con said...

Yes PS has file operators, and you can use them in GhistScript. See http://www.ghostscript.com/doc/current/Use.htm#Finding_files

Saturday, 08 October, 2011

Post a Comment

Subscribe to: Post Comments (Atom)