/test-data/rgtestouts/rgManQQ/rgManQQtest1.html

https://bitbucket.org/cistrome/cistrome-harvard/ · HTML · 255 lines · 217 code · 38 blank · 0 comment · 0 complexity · 40a1212c2e931b506f9235a0beb0aa6c MD5 · raw file 

    

  1. <?xml version="1.0" encoding="utf-8" ?>
    2. 

  2. <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
    3. 

  3. <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
    4. 

  4. <head>
    5. 

  5. <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    6. 

  6. <meta name="generator" content="Galaxy rgManQQ.py tool output - see http://g2.trac.bx.psu.edu/" />
    7. 

  7. <title></title>
    8. 

  8. <link rel="stylesheet" href="/static/style/base.css" type="text/css" />
    9. 

  9. </head>
    10. 

  10. <body>
    11. 

  11. <div class="document">
    12. 

  12. 

  13. <h1>rgManQQtest1</h1>
    14. 

  14. <table>
    15. 

  15. 

  16. <tr><td><a href="Allelep_manhattan.pdf">Allelep_manhattan.pdf</a></td></tr>
    17. 

  17. <tr><td><a href="Allelep_manhattan.pdf"><img src="Allelep_manhattan.png" title="Allelep_manhattan.png" hspace="10" width="800"></a></td></tr>
    18. 

  18. <tr><td><a href="Allelep_qqplot.pdf">Allelep_qqplot.pdf</a></td></tr>
    19. 

  19. <tr><td><a href="Allelep_qqplot.pdf"><img src="Allelep_qqplot.png" title="Allelep_qqplot.png" hspace="10" width="800"></a></td></tr>
    20. 

  20. <tr><td><a href="rgManQQtest1.R">rgManQQtest1.R</a></td></tr>
    21. 

  21. <tr><td><a href="rgManQQtest1.R.log">rgManQQtest1.R.log</a></td></tr>
    22. 

  22. </table>
    23. 

  23. 

  24. <h3>R log follows below</h3><hr><pre>
    25. 

  25. 

  26. Loading required package: reshape
    27. 

  27. 

  28. Loading required package: plyr
    29. 

  29. 

  30. 

  31. 

  32. Attaching package: 'reshape'
    33. 

  33. 

  34. 

  35. 

  36. The following object(s) are masked from 'package:plyr':
    37. 

  37. 

  38. 

  39. 

  40.     rename, round_any
    41. 

  41. 

  42. 

  43. 

  44. Loading required package: grid
    45. 

  45. 

  46. Loading required package: proto
    47. 

  47. 

  48. [1] "### 101 values read from /data/tmp/tmpNaxDwH/database/files/000/dataset_1.dat read - now running plots"
    49. 

  49. 

  50. [1] "## qqplot on Allelep done"
    51. 

  51. 

  52. [1] "## manhattan on Allelep starting 2 3 8"
    53. 

  53. 

  54. [1] "## manhattan plot on Allelep done"
    55. 

  55. 

  56. ## R script=
    57. 

  57. 

  58. # generalised so 3 core fields passed as parameters ross lazarus March 24 2010 for rgenetics
    59. 

  59. # Originally created as qqman with the following 
    60. 

  60. # attribution:
    61. 

  61. #--------------
    62. 

  62. # Stephen Turner
    63. 

  63. # http://StephenTurner.us/
    64. 

  64. # http://GettingGeneticsDone.blogspot.com/
    65. 

  65. 

  66. # Last updated: 19 July 2011 by Ross Lazarus
    67. 

  67. # R code for making manhattan plots and QQ plots from plink output files. 
    68. 

  68. # With GWAS data this can take a lot of memory. Recommended for use on 
    69. 

  69. # 64bit machines only, for now. 
    70. 

  70. 

  71. #
    72. 

  72. 

  73. library(ggplot2)
    74. 

  74. 

  75. coloursTouse = c('firebrick','darkblue','goldenrod','darkgreen')
    76. 

  76. # not too ugly but need a colour expert please...
    77. 

  77. 

  78. 

  79. DrawManhattan = function(pvals=Null,chrom=Null,offset=Null,title=NULL, max.y="max",suggestiveline=0, genomewide=T, size.x.labels=9, 
    80. 

  80.               size.y.labels=10, annotate=F, SNPlist=NULL,grey=0) {
    81. 

  81.         if (annotate & is.null(SNPlist)) stop("You requested annotation but provided no SNPlist!")
    82. 

  82.         genomewideline=NULL # was genomewideline=-log10(5e-8)
    83. 

  83.         n = length(pvals)
    84. 

  84.         if (genomewide) { # use bonferroni since might be only a small region?
    85. 

  85.             genomewideline = -log10(0.05/n) }
    86. 

  86.         offset = as.integer(offset)
    87. 

  87.         if (n > 1000000) { offset = offset/10000 }
    88. 

  88.         else if (n > 10000) { offset = offset/1000}
    89. 

  89.         chro = as.integer(chrom) # already dealt with X and friends?
    90. 

  90.         pvals = as.double(pvals)
    91. 

  91.         d=data.frame(CHR=chro,BP=offset,P=pvals)
    92. 

  92.         if ("CHR" %in% names(d) & "BP" %in% names(d) & "P" %in% names(d) ) {
    93. 

  93.                 d=d[!is.na(d$P), ]
    94. 

  94.                 d=d[!is.na(d$BP), ]
    95. 

  95.                 d=d[!is.na(d$CHR), ]
    96. 

  96.                 #limit to only chrs 1-22, x=23,y=24,Mt=25?
    97. 

  97.                 d=d[d$CHR %in% 1:25, ]
    98. 

  98.                 d=d[d$P>0 & d$P<=1, ]
    99. 

  99.                 d$logp = as.double(-log10(d$P))
    100. 

  100.                 dlen = length(d$P)
    101. 

  101.                 d$pos=NA
    102. 

  102.                 ticks=NULL
    103. 

  103.                 lastbase=0
    104. 

  104.                 chrlist = unique(d$CHR)
    105. 

  105.                 chrlist = as.integer(chrlist)
    106. 

  106.                 chrlist = sort(chrlist) # returns lexical ordering 
    107. 

  107.                 if (max.y=="max") { maxy = ceiling(max(d$logp)) } 
    108. 

  108.                    else { maxy = max.y }
    109. 

  109.                 nchr = length(chrlist) # may be any number?
    110. 

  110.                 maxy = max(maxy,1.1*genomewideline)
    111. 

  111.                 if (nchr >= 2) {
    112. 

  112.                     for (x in c(1:nchr)) {
    113. 

  113.                         i = chrlist[x] # need the chrom number - may not == index
    114. 

  114.                         if (x == 1) { # first time
    115. 

  115.                             d[d$CHR==i, ]$pos = d[d$CHR==i, ]$BP # initialize to first BP of chr1
    116. 

  116.                             dsub = subset(d,CHR==i)
    117. 

  117.                             dlen = length(dsub$P)
    118. 

  118.                             lastbase = max(dsub$pos) # last one
    119. 

  119.                             tks = d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1]
    120. 

  120.                             lastchr = i
    121. 

  121.                         } else {
    122. 

  122.                             d[d$CHR==i, ]$pos = d[d$CHR==i, ]$BP+lastbase # one humongous contig
    123. 

  123.                             if (sum(is.na(lastchr),is.na(lastbase),is.na(d[d$CHR==i, ]$pos))) { 
    124. 

  124.                                 cat(paste('manhattan: For',title,'chrlistx=',i,'lastchr=',lastchr,'lastbase=',lastbase,'pos=',d[d$CHR==i,]$pos))
    125. 

  125.                              }   
    126. 

  126.                             tks=c(tks, d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1])
    127. 

  127.                             lastchr = i
    128. 

  128.                             dsub = subset(d,CHR==i)
    129. 

  129.                             lastbase = max(dsub$pos) # last one
    130. 

  130.                         }
    131. 

  131.                     ticklim=c(min(d$pos),max(d$pos))
    132. 

  132.                     xlabs = chrlist
    133. 

  133.                     }
    134. 

  134.                 } else { # nchr is 1
    135. 

  135.                    nticks = 10
    136. 

  136.                    last = max(d$BP)
    137. 

  137.                    first = min(d$BP)
    138. 

  138.                    tks = c(first)
    139. 

  139.                    t = (last-first)/nticks # units per tick
    140. 

  140.                    for (x in c(1:(nticks))) { 
    141. 

  141.                         tks = c(tks,round(x*t)+first) }
    142. 

  142.                    ticklim = c(first,last)
    143. 

  143.                 } # else
    144. 

  144.                 if (grey) {mycols=rep(c("gray10","gray60"),max(d$CHR))
    145. 

  145.                            } else {
    146. 

  146.                            mycols=rep(coloursTouse,max(d$CHR))
    147. 

  147.                            }
    148. 

  148.                 dlen = length(d$P)
    149. 

  149.                 d$pranks = rank(d$P)/dlen
    150. 

  150.                 d$centiles = 100*d$pranks # small are interesting
    151. 

  151.                 d$sizes = ifelse((d$centile < 1),2,1)
    152. 

  152.                 if (annotate) d.annotate=d[as.numeric(substr(d$SNP,3,100)) %in% SNPlist, ]
    153. 

  153.                 if (nchr >= 2) {
    154. 

  154.                         manplot=qplot(pos,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR),size=factor(sizes))
    155. 

  155.                         manplot=manplot+scale_x_continuous(name="Chromosome", breaks=tks, labels=xlabs,limits=ticklim) 
    156. 

  156.                         manplot=manplot+scale_size_manual(values = c(0.5,1.5)) # requires discreet scale - eg factor
    157. 

  157.                         #manplot=manplot+scale_size(values=c(0.5,2)) # requires continuous 
    158. 

  158.                         }
    159. 

  159.                 else {
    160. 

  160.                         manplot=qplot(BP,logp,data=d, ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
    161. 

  161.                         manplot=manplot+scale_x_continuous(name=paste("Chromosome",chrlist[1]), breaks=tks, labels=tks,limits=ticklim) 
    162. 

  162.                      }                 
    163. 

  163.                 manplot=manplot+scale_y_continuous(limits=c(0,maxy), breaks=1:maxy, labels=1:maxy)
    164. 

  164.                 manplot=manplot+scale_colour_manual(value=mycols)
    165. 

  165.                 if (annotate) {  manplot=manplot + geom_point(data=d.annotate, colour=I("green3")) } 
    166. 

  166.                 manplot=manplot + opts(legend.position = "none") 
    167. 

  167.                 manplot=manplot + opts(title=title)
    168. 

  168.                 manplot=manplot+opts(
    169. 

  169.                      panel.background=theme_blank(), 
    170. 

  170.                      axis.text.x=theme_text(size=size.x.labels, colour="grey50"), 
    171. 

  171.                      axis.text.y=theme_text(size=size.y.labels, colour="grey50"), 
    172. 

  172.                      axis.ticks=theme_segment(colour=NA)
    173. 

  173.                 )
    174. 

  174.                 if (suggestiveline) manplot=manplot+geom_hline(yintercept=suggestiveline,colour="blue", alpha=I(1/3))
    175. 

  175.                 if (genomewideline) manplot=manplot+geom_hline(yintercept=genomewideline,colour="red")
    176. 

  176.                 manplot
    177. 

  177.         }       else {
    178. 

  178.                 stop("Make sure your data frame contains columns CHR, BP, and P")
    179. 

  179.         }
    180. 

  180. }
    181. 

  181. 

  182. 

  183. 

  184. qq = function(pvector, title=NULL, spartan=F) {
    185. 

  185.         # Thanks to Daniel Shriner at NHGRI for providing this code for creating expected and observed values
    186. 

  186.         o = -log10(sort(pvector,decreasing=F))
    187. 

  187.         e = -log10( 1:length(o)/length(o) )
    188. 

  188.         # you could use base graphics
    189. 

  189.         # plot(e,o,pch=19,cex=0.25, xlab=expression(Expected~~-log[10](italic(p))), 
    190. 

  190.         # ylab=expression(Observed~~-log[10](italic(p))), xlim=c(0,max(e)), ylim=c(0,max(e)))
    191. 

  191.         # lines(e,e,col="red")
    192. 

  192.         #You'll need ggplot2 installed to do the rest
    193. 

  193.         qq=qplot(e,o, xlim=c(0,max(e)), ylim=c(0,max(o))) + stat_abline(intercept=0,slope=1, col="red")
    194. 

  194.         qq=qq+opts(title=title)
    195. 

  195.         qq=qq+scale_x_continuous(name=expression(Expected~~-log[10](italic(p))))
    196. 

  196.         qq=qq+scale_y_continuous(name=expression(Observed~~-log[10](italic(p))))
    197. 

  197.         if (spartan) plot=plot+opts(panel.background=theme_rect(col="grey50"), panel.grid.minor=theme_blank())
    198. 

  198.         qq
    199. 

  199. }
    200. 

  200. 

  201. rgqqMan = function(infile="/data/tmp/tmpNaxDwH/database/files/000/dataset_1.dat",chromcolumn=2, offsetcolumn=3, pvalscolumns=c(8), 
    202. 

  202. title="rgManQQtest1",grey=0) {
    203. 

  203. rawd = read.table(infile,head=T,sep='\t')
    204. 

  204. dn = names(rawd)
    205. 

  205. cc = dn[chromcolumn]
    206. 

  206. oc = dn[offsetcolumn] 
    207. 

  207. rawd[,cc] = sub('chr','',rawd[,cc],ignore.case = T) # just in case
    208. 

  208. rawd[,cc] = sub(':','',rawd[,cc],ignore.case = T) # ugh
    209. 

  209. rawd[,cc] = sub('X',23,rawd[,cc],ignore.case = T)
    210. 

  210. rawd[,cc] = sub('Y',24,rawd[,cc],ignore.case = T)
    211. 

  211. rawd[,cc] = sub('Mt',25,rawd[,cc], ignore.case = T)
    212. 

  212. nams = c(cc,oc) # for sorting
    213. 

  213. plen = length(rawd[,1])
    214. 

  214. print(paste('###',plen,'values read from',infile,'read - now running plots',sep=' '))
    215. 

  215. rawd = rawd[do.call(order,rawd[nams]),]
    216. 

  216. # mmmf - suggested by http://onertipaday.blogspot.com/2007/08/sortingordering-dataframe-according.html
    217. 

  217. # in case not yet ordered
    218. 

  218. if (plen > 0) {
    219. 

  219.   for (pvalscolumn in pvalscolumns) {
    220. 

  220.   if (pvalscolumn > 0) 
    221. 

  221.      {
    222. 

  222.      cname = names(rawd)[pvalscolumn]
    223. 

  223.      mytitle = paste('p=',cname,', ',title,sep='')
    224. 

  224.      myfname = chartr(' ','_',cname)
    225. 

  225.      myqqplot = qq(rawd[,pvalscolumn],title=mytitle)
    226. 

  226.      ggsave(filename=paste(myfname,"qqplot.png",sep='_'),myqqplot,width=8,height=6,dpi=96)
    227. 

  227.      ggsave(filename=paste(myfname,"qqplot.pdf",sep='_'),myqqplot,width=8,height=6,dpi=96)
    228. 

  228.      print(paste('## qqplot on',cname,'done'))
    229. 

  229.      if ((chromcolumn > 0) & (offsetcolumn > 0)) {
    230. 

  230.          print(paste('## manhattan on',cname,'starting',chromcolumn,offsetcolumn,pvalscolumn))
    231. 

  231.          mymanplot= DrawManhattan(chrom=rawd[,chromcolumn],offset=rawd[,offsetcolumn],pvals=rawd[,pvalscolumn],title=mytitle,grey=grey)
    232. 

  232.          ggsave(filename=paste(myfname,"manhattan.png",sep='_'),mymanplot,width=8,height=6,dpi=96)
    233. 

  233.          ggsave(filename=paste(myfname,"manhattan.pdf",sep='_'),mymanplot,width=8,height=6,dpi=96)
    234. 

  234.          print(paste('## manhattan plot on',cname,'done'))
    235. 

  235.          }
    236. 

  236.          else {
    237. 

  237.               print(paste('chrom column =',chromcolumn,'offset column = ',offsetcolumn,
    238. 

  238.               'so no Manhattan plot - supply both chromosome and offset as numerics for Manhattan plots if required'))
    239. 

  239.               } 
    240. 

  240.      } 
    241. 

  241.   else {
    242. 

  242.         print(paste('pvalue column =',pvalscolumn,'Cannot parse it so no plots possible'))
    243. 

  243.       }
    244. 

  244.   } # for pvalscolumn
    245. 

  245.  } else { print('## Problem - no values available to plot - was there really a chromosome and offset column?') }
    246. 

  246. }
    247. 

  247. 

  248. rgqqMan() 
    249. 

  249. # execute with defaults as substituted
    250. 

  250. 

  251. </pre>
    252. 

  252. 

  253. <h3><a href="http://rgenetics.org">Rgenetics</a> tool rgManQQ.py run at 04/08/2011 11:46:22</h3>
    254. 

  254. </div></body></html>
    255.