2016년 5월 11일 수요일

gnuplot으로 위 아래가 붙은 그림 그리기

X 축의 데이터가 동일한 두 장의 플롯을 여백 없이 세로로 이어붙이면 큰 면적을 차지하지 않으면서 효과적인 그림을 나타낼 수 있다. 전문적인 용어를 쓰자면 multiplot의 응용 사례중 하나로서 "stacked plots"이 되겠다. 다음 링크를 연구하면 해결 방법이 나온다.

gnuplot demo script: layout.dem (gnuplot 5.0)

CentOS 배포판에서 설치되는 gnuplot은 버전이 너무 구식이다. CentOS 6.7은 아직도 version 4.2 patchlevel 6이다. 이를 version 5.0으로 업그레이드하니 그림도 더욱 예뻐지고 GUI 형식의 조절 창도 열린다. 문제는 업그레이드를 하고 나니 MUMmer에서 만든 스크립트에서 다음의 명령을 이해하지 못한다는 것이다.

set mouse clipboardformat "[%.0f, %.0f]"
          ^
"T.gp", line 94: wrong option
버전 5.0의 매뉴얼에서 set mouse 명령의 유효한 옵션이 어떻게 바뀌었는지 차차 알아보기로 하자. 가장 단순무식한 해결 방법은 gnuplot script에서 해당 라인을 지워버리면 그만이다. 어차피 마우스를 인터랙티브하게 사용할 것은 아니므로.

오늘 그리려는 그림은 single reference genome 서열 위에 query sequence(multiple)을 MUMmer로 정렬하여 순서를 맞춘 그림 2장을 위 아래로 붙여서 표현한 것이다. 다시 말하자면 하나의 reference sequence에 대해서 두 종류의 contig set을 각각 정렬한 결과를 비교하자는 뜻이다. Stacked plots의 기본 트릭은 총 3장의 plot이 필요하다면 4, 1 레이아웃으로 선언을 하되 plot 사이의 마진은 없애고, 4번째 plot이 위치할 곳에 X 축과 설명을 붙이는 것이다. 문제는 query sequence가 많으면 Y축 라벨이 지저분하게 나타나므로 이를 취향에 맞도록 적절히 조절하는 것이다. 이번의 경우에는 ytic을 숫자만 없이 표시하게 하였다. 이것 외에도 몇가지의 트릭이 들어갔다.

  • Y 라벨에서 밑줄(_)을 제대로 표현하기: 이중으로 이스케이프를 하지 않으면 아래첨자가 된다.
  • 특정 위치에 수직선 세우기: reference 서열은 실제로는 염색체 + 플라스미드를 하나로 연결한 것이다. 경계가 되는 곳에 빨강색 세로 화살표(headless)를 세웠다.
  • X축의 표시 단위를 지수형태로 만들었다. (format specifier 참조)
단일 reference에 대하여 multiple query sequence를 정렬하는 스크립트로는 run-mummer4.sh라는 것이 있다. 어디서 났느냐고? 내가 만들었다.


#!/bin/sh
#
# run-mummer4.sh: MUMmer-based shell script for aligning draft sequences #                        to a finished genome.
#
if [ $# = 0 ]
    then
        echo "Usage : run-mummer4.sh [-q]" 1>&2
        echo "The last option is for delta-filter."
        echo "Other nucmer default options: -maxmatch -c 100"
          exit 1
fi
#nucmer --prefix=$1 -mum -c 100 $2 $3
nucmer --prefix=$1 -mum -c 1000 $2 $3
if [ "$4" = "-q" ]
    then
        delta-filter $4 $1.delta > tmp
        mv $1.delta $1.delta.org
        mv tmp $1.delta
fi
show-coords -rcl $1.delta > $1.coords
REFNAME=`awk '$1~/^>/{print $1}' $2 | sed 's/^>//'`
for QRYNAME in `awk '$1~/^>/{print $1}' $3 | sed 's/^>//'`
do
    show-aligns $1.delta $REFNAME $QRYNAME > $1.aligns
done
mummerplot --postscript $1.delta -R $2 -Q $3  --layout --prefix=$1
이를 (reference + query 1), (refernece + query 2)로 각각 실행하여 gnuplot 스크립트를 만든 뒤 이어붙여서 다음과 같이 편집하였다. 위에서 소개한 트릭들이 어떻게 반영되었는지를 잘 확인해 보라. ytics를 선언한 곳은 너무 길어서 적당히 생략하였다.

set multiplot layout 3,1 title "Auto-layout of stacked plots\n" font ",12"
set grid
unset title
unset key
set border
set ylabel "DS1"
set tmargin 0
set bmargin 0
set lmargin 10
set rmargin 3
set xtics format " "
set ytics ( \
 "" 1, \
 "" 36204, \
...
 "" 4969774, \
 "" 4971242 \
)
#set tics scale 1
set xrange [1:5848607]
set yrange [1:4971242]
set style line 1  lt 1 lw 2 pt 6 ps 0.5
set style line 2  lt 3 lw 2 pt 6 ps 0.5
set style line 3  lt 2 lw 2 pt 6 ps 0.5
set arrow from 5461741,0 to 5461741,4971242 nohead lc "red"
plot \
 "DS1.fplot" title "FWD" w lp ls 1, \
 "DS1.rplot" title "REV" w lp ls 2
set ytics ( \
 "" 1, \
 "" 17170, \
...
 "" 5880317, \
 "" 5882187 \
)
set xlabel "Chromosome (bp)"
set ylabel "2\\_A\\_57\\_GC2"
set format x "%.0tx10^%1T"
set xtics
set xrange [1:5848607]
set yrange [1:5882187]
set style line 1  lt 1 lw 2 pt 6 ps 0.5
set style line 2  lt 3 lw 2 pt 6 ps 0.5
set style line 3  lt 2 lw 2 pt 6 ps 0.5
set arrow from 5461741,0 to 5461741,5882178 nohead lc "red"
plot \
 "2_A.fplot" title "FWD" w lp ls 1, \
 "2_A.rplot" title "REV" w lp ls 2
pause -1
결과물은 다음과 같다.







댓글 없음: