root/galaxy-central/test-data/rgtestouts/rgManQQ/rgManQQtest1.html @ 2

リビジョン 2, 9.8 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central
行番号 
1<?xml version="1.0" encoding="utf-8" ?>
2<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
3<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
4<head>
5<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
6<meta name="generator" content="Galaxy rgManQQ.py tool output - see http://g2.trac.bx.psu.edu/" />
7<title></title>
8<link rel="stylesheet" href="/static/style/base.css" type="text/css" />
9</head>
10<body>
11<div class="document">
12
13<h1>rgManQQtest1</h1>
14<table>
15
16<tr><td><a href="Allelep_manhattan.png"><img src="Allelep_manhattan.png" alt="Allelep_manhattan.png hspace="10" width="400"><br>(Click to download image Allelep_manhattan.png)</a></td></tr>
17<tr><td><a href="Allelep_qqplot.png"><img src="Allelep_qqplot.png" alt="Allelep_qqplot.png hspace="10" width="400"><br>(Click to download image Allelep_qqplot.png)</a></td></tr>
18<tr><td><a href="rgManQQtest1.R">rgManQQtest1.R</a></td></tr>
19<tr><td><a href="rgManQQtest1.R.log">rgManQQtest1.R.log</a></td></tr>
20</table>
21
22<h3>R log follows below</h3><hr><pre>
23
24Loading required package: reshape
25
26Loading required package: plyr
27
28
29
30Attaching package: 'reshape'
31
32
33
34The following object(s) are masked from 'package:plyr':
35
36
37
38    round_any
39
40
41
42Loading required package: grid
43
44Loading required package: proto
45
46[1] "### 101 values read from /tmp/rgManQQtempcWfFkc read - now running plots"
47
48[1] "## qqplot on Allelep done"
49
50[1] "## manhattan on Allelep starting 1 2 3"
51
52[1] "## manhattan plot on Allelep done"
53
54## R script=
55
56# license not stated so I'm assuming LGPL is ok for my derived work?
57# generalised so 3 core fields passed as parameters ross lazarus March 24 2010 for rgenetics
58# Originally created as qqman with the following
59# attribution:
60#--------------
61# Stephen Turner
62# http://StephenTurner.us/
63# http://GettingGeneticsDone.blogspot.com/
64
65# Last updated: Tuesday, December 22, 2009
66# R code for making manhattan plots and QQ plots from plink output files.
67# With GWAS data this can take a lot of memory. Recommended for use on
68# 64bit machines only, for now.
69
70#
71
72library(ggplot2)
73
74coloursTouse = c('firebrick','darkblue','goldenrod','darkgreen')
75# not too fugly but need a colour expert please...
76
77
78manhattan = function(chrom=NULL,offset=NULL,pvals=NULL, title=NULL, max.y="max",
79   suggestiveline=0, genomewide=T, size.x.labels=9, size.y.labels=10, annotate=F, SNPlist=NULL,grey=0) {
80
81        if (annotate & is.null(SNPlist)) stop("You requested annotation but provided no SNPlist!")
82        genomewideline=NULL # was genomewideline=-log10(5e-8)
83        if (genomewide) { # use bonferroni since might be only a small region?
84            genomewideline = -log10(0.05/length(pvals)) }
85        d=data.frame(CHR=chrom,BP=offset,P=pvals)
86
87        #limit to only chrs 1-23?
88        d=d[d$CHR %in% 1:23, ]
89
90        if ("CHR" %in% names(d) & "BP" %in% names(d) & "P" %in% names(d) ) {
91                d=na.omit(d)
92                d=d[d$P>0 & d$P<=1, ]
93                d$logp = -log10(d$P)
94
95                d$pos=NA
96                ticks=NULL
97                lastbase=0
98                chrlist = unique(d$CHR)
99                nchr = length(chrlist) # may be any number?
100                if (nchr >= 2) {
101                for (x in c(1:nchr)) {
102                        i = chrlist[x] # need the chrom number - may not == index
103                        if (x == 1) { # first time
104                                d[d$CHR==i, ]$pos=d[d$CHR==i, ]$BP
105                                tks = d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1]
106                        }       else {
107                                lastchr = chrlist[x-1] # previous whatever the list
108                                lastbase=lastbase+tail(subset(d,CHR==lastchr)$BP, 1)
109                                d[d$CHR==i, ]$pos=d[d$CHR==i, ]$BP+lastbase
110                                tks=c(tks, d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1])
111                        }
112                    ticklim=c(min(d$pos),max(d$pos))
113                    xlabs = chrlist
114                    }
115                } else { # nchr is 1
116                   nticks = 10
117                   last = max(offset)
118                   first = min(offset)
119                   tks = c()
120                   t = (last-first)/nticks # units per tick
121                   for (x in c(1:nticks)) tks = c(tks,round(x*t))
122                   xlabs = tks
123                   ticklim = c(first,last)
124                } # else
125                if (grey) {mycols=rep(c("gray10","gray60"),max(d$CHR))
126                           } else {
127                           mycols=rep(coloursTouse,max(d$CHR))
128                           }
129
130                if (max.y=="max") maxy=ceiling(max(d$logp)) else maxy=max.y
131                maxy = max(maxy,1.1*genomewideline)
132                # if (maxy<8) maxy=8
133                # only makes sense if genome wide is assumed - we could have a fine mapping region? 
134                if (annotate) d.annotate=d[as.numeric(substr(d$SNP,3,100)) %in% SNPlist, ]
135                if (nchr >= 2) {
136                        manplot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
137                        manplot=manplot+scale_x_continuous(name="Chromosome", breaks=tks, labels=xlabs) }
138                else {
139                        manplot=qplot(BP,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
140                        manplot=manplot+scale_x_continuous("BP") }                 
141                manplot=manplot+scale_y_continuous(limits=c(0,maxy), breaks=1:maxy, labels=1:maxy)
142                manplot=manplot+scale_colour_manual(value=mycols)
143                if (annotate) {  manplot=manplot + geom_point(data=d.annotate, colour=I("green3")) }
144                manplot=manplot + opts(legend.position = "none")
145                manplot=manplot + opts(title=title)
146                manplot=manplot+opts(
147                        panel.background=theme_blank(),
148                        axis.text.x=theme_text(size=size.x.labels, colour="grey50"),
149                        axis.text.y=theme_text(size=size.y.labels, colour="grey50"),
150                        axis.ticks=theme_segment(colour=NA)
151                )
152                #manplot = manplot + opts(panel.grid.y.minor=theme_blank(),panel.grid.y.major=theme_blank())
153                #manplot = manplot + opts(panel.grid.major=theme_blank())
154                 
155                if (suggestiveline) manplot=manplot+geom_hline(yintercept=suggestiveline,colour="blue", alpha=I(1/3))
156                if (genomewideline) manplot=manplot+geom_hline(yintercept=genomewideline,colour="red")
157                manplot
158        }       else {
159                stop("Make sure your data frame contains columns CHR, BP, and P")
160        }
161}
162
163
164
165qq = function(pvector, title=NULL, spartan=F) {
166        # Thanks to Daniel Shriner at NHGRI for providing this code for creating expected and observed values
167        o = -log10(sort(pvector,decreasing=F))
168        e = -log10( 1:length(o)/length(o) )
169        # you could use base graphics
170        # plot(e,o,pch=19,cex=0.25, xlab=expression(Expected~~-log[10](italic(p))),
171        # ylab=expression(Observed~~-log[10](italic(p))), xlim=c(0,max(e)), ylim=c(0,max(e)))
172        # lines(e,e,col="red")
173        #You'll need ggplot2 installed to do the rest
174        qq=qplot(e,o, xlim=c(0,max(e)), ylim=c(0,max(o))) + stat_abline(intercept=0,slope=1, col="red")
175        qq=qq+opts(title=title)
176        qq=qq+scale_x_continuous(name=expression(Expected~~-log[10](italic(p))))
177        qq=qq+scale_y_continuous(name=expression(Observed~~-log[10](italic(p))))
178        if (spartan) plot=plot+opts(panel.background=theme_rect(col="grey50"), panel.grid.minor=theme_blank())
179        qq
180}
181rgqqMan = function(infile="/tmp/rgManQQtempcWfFkc",chromcolumn=1, offsetcolumn=2, pvalscolumns=c(3),
182title="rgManQQtest1",grey=0) {
183rawd = read.table(infile,head=T,sep='\t')
184dn = names(rawd)
185cc = dn[chromcolumn]
186oc = dn[offsetcolumn]
187nams = c(cc,oc)
188plen = length(rawd[,1])
189doreorder=1
190print(paste('###',plen,'values read from',infile,'read - now running plots',sep=' '))
191if (plen > 0) {
192  for (pvalscolumn in pvalscolumns) {
193  if (pvalscolumn > 0)
194     {
195     cname = names(rawd)[pvalscolumn]
196     mytitle = paste('p=',cname,', ',title,sep='')
197     myfname = chartr(' ','_',cname)
198     myqqplot = qq(rawd[,pvalscolumn],title=mytitle)
199     ggsave(filename=paste(myfname,"qqplot.png",sep='_'),myqqplot,width=11,height=8,dpi=100)
200     print(paste('## qqplot on',cname,'done'))
201     if ((chromcolumn > 0) & (offsetcolumn > 0)) {
202         if (doreorder) {
203             rawd = rawd[do.call(order,rawd[nams]),]
204             # mmmf - suggested by http://onertipaday.blogspot.com/2007/08/sortingordering-dataframe-according.html
205             # in case not yet ordered
206             doreorder = 0
207             }
208         print(paste('## manhattan on',cname,'starting',chromcolumn,offsetcolumn,pvalscolumn))
209         mymanplot= manhattan(chrom=rawd[,chromcolumn],offset=rawd[,offsetcolumn],pvals=rawd[,pvalscolumn],title=mytitle,grey=grey)
210         print(paste('## manhattan plot on',cname,'done'))
211         ggsave(filename=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=11,height=8,dpi=100)
212         }
213         else {
214              print(paste('chrom column =',chromcolumn,'offset column = ',offsetcolumn,
215              'so no Manhattan plot - supply both chromosome and offset as numerics for Manhattan plots if required'))
216              }
217     }
218  else {
219        print(paste('pvalue column =',pvalscolumn,'Cannot parse it so no plots possible'))
220      }
221  } # for pvalscolumn
222 } else { print('## Problem - no values available to plot - was there really a chromosome and offset column?') }
223}
224
225rgqqMan()
226# execute with defaults as substituted
227
228</pre>
229
230<h3><a href="http://rgenetics.org">Rgenetics</a> tool rgManQQ.py run at 18/09/2010 20:48:26</h3>
231</div></body></html>
232
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。