Mercurial > repos > mvdbeek > r_goseq_1_22_0

comparison get_length_and_gc_content.r @ 1:3ab168143b69 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 05:26:51 -0500
parents
children f9b964d1d386
comparison
equal deleted inserted replaced
0:fe71b97cc1a5 1:3ab168143b69
1 # originally by Devon Ryan, https://www.biostars.org/p/84467/
2
3 library(GenomicRanges)
4 library(rtracklayer)
5 library(Rsamtools)
6 library(optparse)
7
8 option_list <- list(
9 make_option(c("-g","--gtf"), type="character", help="Input GTF file with gene / exon information."),
10 make_option(c("-f","--fasta"), type="character", default=FALSE, help="Fasta file that corresponds to the supplied GTF."),
11 make_option(c("-o","--output"), type="character", default=FALSE, help="Output file with gene name, length and GC content.")
12 )
13
14 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
15 args = parse_args(parser)
16
17 GTFfile = args.gtf
18 FASTAfile = args.fasta
19 output = args.output
20
21 #Load the annotation and reduce it
22 GTF <- import.gff(GTFfile, format="gtf", genome=NA, asRangedData=F, feature.type="exon")
23 grl <- reduce(split(GTF, elementMetadata(GTF)$gene_name))
24 reducedGTF <- unlist(grl, use.names=T)
25 elementMetadata(reducedGTF)$gene_name <- rep(names(grl), elementLengths(grl))
26
27 #Open the fasta file
28 FASTA <- FaFile(FASTAfile)
29 open(FASTA)
30
31 #Add the GC numbers
32 elementMetadata(reducedGTF)$nGCs <- letterFrequency(getSeq(FASTA, reducedGTF), "GC")[,1]
33 elementMetadata(reducedGTF)$widths <- width(reducedGTF)
34
35 #Create a list of the ensembl_id/GC/length
36 calc_GC_length <- function(x) {
37 nGCs = sum(elementMetadata(x)$nGCs)
38 width = sum(elementMetadata(x)$widths)
39 c(width, nGCs/width)
40 }
41 output <- t(sapply(split(reducedGTF, elementMetadata(reducedGTF)$gene_name), calc_GC_length))
42 colnames(output) <- c("Length", "GC")
43
44 write.table(output, file="GC_lengths.tsv", sep="\t")