Mercurial > repos > mvdbeek > r_goseq_1_22_0
comparison get_length_and_gc_content.r @ 1:3ab168143b69 draft
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") |