Mercurial > repos > mvdbeek > r_goseq_1_22_0

comparison goseq.r @ 16:8ce951313688 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 Thu, 25 Feb 2016 08:49:20 -0500
parents fe71b97cc1a5
children 1b03f6232900
comparison
equal deleted inserted replaced
15:e3b1dd2f7f70 16:8ce951313688
11 make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=FALSE, help="Path to plot comparing sampling with wallenius p-values."), 11 make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=FALSE, help="Path to plot comparing sampling with wallenius p-values."),
12 make_option(c("-padj", "--p_adj_column"), type="integer",help="Column that contains p. adjust values"), 12 make_option(c("-padj", "--p_adj_column"), type="integer",help="Column that contains p. adjust values"),
13 make_option(c("-c", "--cutoff"), type="double",dest="p_adj_cutoff", 13 make_option(c("-c", "--cutoff"), type="double",dest="p_adj_cutoff",
14 help="Genes with p.adjust below cutoff are considered not differentially expressed and serve as control genes"), 14 help="Genes with p.adjust below cutoff are considered not differentially expressed and serve as control genes"),
15 make_option(c("-r", "--repcnt"), type="integer", default=100, help="Number of repeats for sampling"), 15 make_option(c("-r", "--repcnt"), type="integer", default=100, help="Number of repeats for sampling"),
16 make_option(c("-lf", "--length_file"), default=FALSE, help = "Path to "),
16 make_option(c("-g", "--genome"), type="character", help = "Genome [used for looking up correct gene length]"), 17 make_option(c("-g", "--genome"), type="character", help = "Genome [used for looking up correct gene length]"),
17 make_option(c("-i", "--gene_id"), type="character", help="Gene ID of gene column in DGE file") 18 make_option(c("-i", "--gene_id"), type="character", help="Gene ID of gene column in DGE file")
18 ) 19 )
19 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list) 20 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
20 args = parse_args(parser) 21 args = parse_args(parser)
21 22
22 # Vars: 23 # Vars:
23 dge_file = args$dge_file 24 dge_file = args$dge_file
24 p_adj_column = args$p_adj_colum 25 p_adj_column = args$p_adj_colum
25 p_adj_cutoff = args$p_adj_cutoff 26 p_adj_cutoff = args$p_adj_cutoff
27 length_file = args$length_file
26 genome = args$genome 28 genome = args$genome
27 gene_id = args$gene_id 29 gene_id = args$gene_id
28 wallenius_tab = args$wallenius_tab 30 wallenius_tab = args$wallenius_tab
29 sampling_tab = args$sampling_tab 31 sampling_tab = args$sampling_tab
30 nobias_tab = args$nobias_tab 32 nobias_tab = args$nobias_tab
36 # format DE genes into vector suitable for use with goseq 38 # format DE genes into vector suitable for use with goseq
37 dge_table = read.delim(dge_file, header = TRUE, sep="\t", check.names = FALSE) 39 dge_table = read.delim(dge_file, header = TRUE, sep="\t", check.names = FALSE)
38 genes = as.integer(dge_table[,p_adj_column]<p_adj_cutoff) 40 genes = as.integer(dge_table[,p_adj_column]<p_adj_cutoff)
39 names(genes) = dge_table[,1] # Assuming first row contains gene names 41 names(genes) = dge_table[,1] # Assuming first row contains gene names
40 42
43 # Get gene lengths
44
45 if (length_file) {
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
48 } else {
49 gene_lengths = getlength(names(genes),genome, gene_id)
50 }
51
41 # Estimate PWF 52 # Estimate PWF
42 53
54 print(gene_lengths)
55
43 pdf(length_bias_plot) 56 pdf(length_bias_plot)
44 pwf=nullp(genes, genome , gene_id) 57 pwf=nullp(genes, genome, gene_id, gene_lengths)
45 dev.off() 58 dev.off()
46 # Null dstribution wallenius 59 # Null dstribution wallenius
47 GO.wall=goseq(pwf, genome, gene_id) 60 GO.wall=goseq(pwf, genome, gene_id)
48 61
49 GO.nobias=goseq(pwf, genome, gene_id, method="Hypergeometric") 62 GO.nobias=goseq(pwf, genome, gene_id, method="Hypergeometric")
62 75
63 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)
64 write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE) 77 write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
65 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE) 78 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
66 79
80 sessionInfo()
81
67 # Use the following to get a list of supported genomes / gene ids 82 # Use the following to get a list of supported genomes / gene ids
68 83
69 # write.table(supportedGenomes(), "available_genomes.tab", row.names = FALSE, quote=FALSE) 84 # write.table(supportedGenomes(), "available_genomes.tab", row.names = FALSE, quote=FALSE)
70 # write.table(supportedGeneIDs(), "supported_gene_ids.tab", row.name = FALSE, quote = FALSE) 85 # write.table(supportedGeneIDs(), "supported_gene_ids.tab", row.name = FALSE, quote = FALSE)
71 # write.table(table.summary, "input_gene_count_matrix.tab", row.names = FALSE, quote = FALSE) 86 # write.table(table.summary, "input_gene_count_matrix.tab", row.names = FALSE, quote = FALSE)