Mercurial > repos > mvdbeek > r_goseq_1_22_0

comparison goseq.r @ 17:1b03f6232900 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/goseq_1_22_0 commit fdd0811efc61c31f88ff17096fbe8ee8cfacd766-dirty
author mvdbeek
date Fri, 26 Feb 2016 08:49:59 -0500
parents 8ce951313688
children 5fb82111ec62
comparison
equal deleted inserted replaced
16:8ce951313688 17:1b03f6232900
44 44
45 if (length_file) { 45 if (length_file) {
46 length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE, row.names=TRUE) 46 length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE, row.names=TRUE)
47 gene_lengths = length_table[names(genes),]$length 47 gene_lengths = length_table[names(genes),]$length
48 } else { 48 } else {
49 gene_lengths = getlength(names(genes),genome, gene_id) 49 gene_lengths = getlength(names(genes), genome, gene_id)
50 } 50 }
51 51
52 # Estimate PWF 52 # Estimate PWF
53 53
54 print(gene_lengths)
55
56 pdf(length_bias_plot) 54 pdf(length_bias_plot)
57 pwf=nullp(genes, genome, gene_id, gene_lengths) 55 pwf=nullp(genes, genome, gene_id, gene_lengths)
58 dev.off() 56 dev.off()
59 # Null dstribution wallenius 57 # wallenius approximation of p-values
60 GO.wall=goseq(pwf, genome, gene_id) 58 GO.wall=goseq(pwf, genome, gene_id)
61 59
62 GO.nobias=goseq(pwf, genome, gene_id, method="Hypergeometric") 60 GO.nobias=goseq(pwf, genome, gene_id, method="Hypergeometric")
63 61
64 # Sampling dsitribution 62 # Sampling distribution
65 GO.samp=goseq(pwf,genome, gene_id, method="Sampling",repcnt=repcnt) 63 if (repcnt > 0) {
66 64 GO.samp=goseq(pwf,genome, gene_id, method="Sampling", repcnt=repcnt)
67 # Compare sampling with wallenius 65 # Compare sampling with wallenius
68 pdf(sample_vs_wallenius_plot) 66 pdf(sample_vs_wallenius_plot)
69 plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]), 67 plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]),
70 xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)", 68 xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)",
71 xlim=c(-3,0)) 69 xlim=c(-3,0))
72 abline(0,1,col=3,lty=2) 70 abline(0,1,col=3,lty=2)
73 dev.off() 71 dev.off()
72 write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
73 }
74 74
75 75
76 write.table(GO.wall, wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE) 76 write.table(GO.wall, wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
77 write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
78 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE) 77 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
79 78
80 sessionInfo() 79 sessionInfo()
81 80
82 # Use the following to get a list of supported genomes / gene ids 81 # Use the following to get a list of supported genomes / gene ids