| Previous changeset 25:ea3f7a6504e0 (2014-12-10) |
|
Commit message:
Uploaded |
|
added:
RNAseqDataAnnotation/RNAseqDataAnnotation.R RNAseqDataAnnotation/RNAseqDataAnnotation.xml |
| b |
| diff -r ea3f7a6504e0 -r f183f8648c5a RNAseqDataAnnotation/RNAseqDataAnnotation.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RNAseqDataAnnotation/RNAseqDataAnnotation.R Wed Dec 10 06:42:21 2014 -0500 |
| [ |
| b'@@ -0,0 +1,223 @@\n+#Author : keime / lornage\n+#Date : 2014/11\n+\n+\n+########################################################################################################\n+#This function concatenates htseq-count result files, normalizes data and annotates data using Ensembl annotations\n+\n+#arguments\n+#Species : Name of the species\n+#ensversion : version of Ensembl to use\n+#fileout : tab-delimited file containing for each library ; gene id, raw read counts, normalized data as well as normalized data/gene length\n+#corresp : data.frame linking the file loaded into galaxy to the corresponding condition\n+#nfiles : number of files(conditions)\n+\n+#output : a data.frame with the following columns :\n+#ensembl gene id\n+#raw read counts for each library (one column per library)\n+#normalized data for each library (one column per library) \n+#normalized data divided by gene length for each library (one column per library)\n+#Gene name\n+#Description\n+\n+#require : biomaRt and DESeq2 Bioconductor packages / package plyr1.8.1\n+\n+#Methods : \n+#Normalization is performed using the method described in Genome Biology 2010;11(10):R106 \n+#and implemented in the DESeq2 Bioconductor package\n+#Gene length correspond to the median of the size of all transcripts corresponding to this gene\n+#########################################################################################################\n+\n+\n+\n+RNAseqDataAnnotation = function(Species, ensversion, fileout, corresp ,nfiles){\t\t\t\t\t\t\t\t\t\t\t\n+\t\n+ #Create a string vector called libnames that contains the name of the samples \n+\tlibnames=rep("",nfiles)\n+\tfor (i in 1:nfiles){\n+\t\tlibnames[i]=toString(corresp$Sample_name[i])\n+\t} \n+\n+ #For all files in corresp read the corresponding file into R\n+\tsuppressPackageStartupMessages(library(plyr, lib.loc = NULL, character.only = FALSE, logical.return = FALSE, warn.conflicts = FALSE, verbose=FALSE, quietly = TRUE))\n+\tdatalist = list()\n+\tfor(i in 1:nfiles){\n+\t\trawdata=read.table(toString(corresp$Files[i]), header =T, sep ="\\t")\n+\t\t#noSpikes_htseq.\n+\t\tnbrrows=nrow(rawdata)\n+\t\tdatalist[[i]]=rawdata[1:(nbrrows-5), ] # skip the last 5 lines of HTSeq-count files\n+\t\tcolnames(datalist[[i]]) = c("ID",libnames[i])\t\t\n+\t} \n+\t\t\n+ #Join all the files in a data.frame called datafile with rownames = gene id\n+\tdatafile = join_all(datalist, by = "ID", type = "left", match = "all")\n+\t\n+ #Calculate the number of geneID pro file\n+\tnbID=data.frame(rep("",nfiles))\n+\tfor(i in 1:nfiles){\n+\t\tnbID[,i]=nrow(datalist[[i]])\n+\t}\n+\ttotalnbID=apply((nbID[,1:nfiles]),1,sum)\n+\t\n+ #Verify that all the files contain the same gene ID\n+\tif (nrow(datafile)*nfiles==totalnbID[1]){\n+ \n+ #Suppress genes not expressed in all samples \n+\t\tdatafile = datafile[apply(datafile[,2:(nfiles+1)],1,sum)!=0,]\n+\t\trow.names(datafile)=datafile[,1]\n+\t\tdata=datafile[,-1]\n+\n+ #Number of libraries\n+\t\tnblib= dim(data)[2]\t\n+ #Determine Data + normalization if the specie is not known \n+\t\tif (Species=="None"){\n+ #Normalized data calculation\n+\t\t\tnbcol = dim(data)[2] #nb of column in the data.frame\n+\t\t\tsuppressPackageStartupMessages(library(DESeq2, lib.loc = NULL, character.only = FALSE, logical.return = FALSE, warn.conflicts = FALSE, verbose=FALSE, quietly = TRUE))\n+\t\t\tconds = factor(1:nblib)\n+\t\t\tdesign = data.frame(Condition=conds)\n+\t\t\tdds = DESeqDataSetFromMatrix(countData=data, colData=design, design=~Condition)\n+\t\t\tdds = estimateSizeFactors(dds)\n+\t\t\tdatanorm = t(t(data)/sizeFactors(dds))\n+\t\t\t\n+ #Data + normalization \n+\t\t\tdataall = data.frame(row.names(datafile), data, datanorm )\n+\t\n+ #Renames columns\n+\t\t\tcolnames(dataall) = c("Ensembl gene id", paste(libnames,"(raw read counts)"), paste(libnames,"(normalized)"))\n+\t\t\twrite.table(dataall, file=fileout, sep="\\t", quote=F, row.names=F)\n+\t\t}\n+ \n+ #Determine Data + normalization + annotation if the specie is known \n+\t\tel'..b'second step \n+\t\t\tnot = rownames(data)[!rownames(data) %in% unique(annotation2$ensembl_gene_id)]\n+\t\t\tif (length(not) !=0){\n+\t\t\t\tif (toString(ensversion)=="oct2014" | toString(ensversion)=="aug2014" ) {\n+\t\t\t\t\tannotationnot = getBM(attributes=c("ensembl_gene_id","external_gene_name","description", "ensembl_transcript_id","exon_chrom_start","exon_chrom_end"),filters="ensembl_gene_id", values=not, mart=ensembl)\n+\t\t\t\t\tannotation = rbind(annotation2, annotationnot)\n+\t\t\t\t} \n+\t\t\t\telse {\n+\t\t\t\t\tannotationnot = getBM(attributes=c("ensembl_gene_id","external_gene_id","description", "ensembl_transcript_id","exon_chrom_start","exon_chrom_end"), filters="ensembl_gene_id", values=not, mart=ensembl)\n+\t\t\t\t\tannotation = rbind(annotation2, annotationnot)\t\t\n+\t\t\t\t}\n+\t\t\t}\t\n+\t\t\telse{\n+\t\t\t\tannotation = annotation2\n+\t\t\t}\n+\t\n+ #Exon length\n+\t\t\tensinfos.exlen = data.frame(annotation$ensembl_gene_id, annotation$ensembl_transcript_id, abs(annotation$exon_chrom_start - annotation$exon_chrom_end)+1)\n+\t\t\tcolnames(ensinfos.exlen) = c("ensembl_gene_id", "ensembl_transcript_id", "exon_length")\n+\t\n+ #Transcript length\n+\t\t\ttlen = tapply(ensinfos.exlen$exon_length, ensinfos.exlen$ensembl_transcript_id, sum)\n+\t\t\ttlen.gene = merge(tlen, unique(ensinfos.exlen[,1:2]), by.x="row.names", by.y="ensembl_transcript_id")\n+\t\t\tcolnames(tlen.gene) = c("ensembl_transcript_id", "transcript_length","ensembl_gene_id")\n+\t\n+ #Gene length = median of the size of all transcripts corresponding to this gene\n+\t\t\tglen = tapply(tlen.gene$transcript_length, tlen.gene$ensembl_gene_id, median)\n+\t\n+ #Data with gene length\n+\t\t\tdatalen = merge(data, glen, by="row.names") \n+\t\t\tcolnames(datalen) = c("Ensembl_gene_id",colnames(data), "Gene_length")\n+\t\n+ #Data with annotations and gene length\n+\t\t\tannotationgene = unique(annotation[,1:3])\n+\t\t\tdataannot = merge(datalen, annotationgene, by.x="Ensembl_gene_id", by.y="ensembl_gene_id")\n+\t\n+ #To keep only the first part of the gene description (before [)\n+\t\t\ttmpdesc = strsplit(as.character(dataannot$description),"[", fixed=T)\n+\t\t\tf = function(l){\n+\t\t\t\tif (length(l)>=1){\n+\t\t\t\t\treturn(l[[1]])\n+\t\t\t\t}\n+\t\t\t\telse{\n+\t\t\t\t\treturn("")\n+\t\t\t\t}\n+\t\t\t}\n+\t\t\ttmpdescok = unlist(lapply(tmpdesc, f))\n+\t\t\tdataannot$description = tmpdescok\n+\t\n+ #Normalized data calculation\n+\t\t\tnbcol = dim(dataannot)[2] #nb of column in the data.frame\n+\t\t\tsuppressPackageStartupMessages(library(DESeq2, lib.loc = NULL, character.only = FALSE, logical.return = FALSE, warn.conflicts = FALSE,verbose=FALSE, quietly = TRUE))\n+\t\t\tconds = factor(1:nblib)\n+\t\t\tdesign = data.frame(Condition=conds)\n+\t\t\tdds = DESeqDataSetFromMatrix(countData=dataannot[,-c(1,nbcol,nbcol-1,nbcol-2)], colData=design, design=~Condition)\n+\t\t\tdds = estimateSizeFactors(dds)\n+\t\t\tdatanorm = t(t(dataannot[,-c(1,nbcol,nbcol-1,nbcol-2)])/sizeFactors(dds))\n+\t\n+ #Normalized data adjusted for gene length (normalized data / gene length)\n+\t\t\trpkn = datanorm / (as.vector(dataannot[,nbcol-2]/1000 ))\n+\t\n+ #Data + annotations + rpkn\n+\t\t\tdataall = data.frame(dataannot[,-c(nbcol,nbcol-1,nbcol-2)] , datanorm, rpkn, dataannot[,c(nbcol-1,nbcol)] )\n+\t\t\n+ #Renames columns\n+\t\t\tcolnames(dataall) = c("Ensembl gene id", paste(libnames,"(raw read counts)"), paste(libnames,"(normalized)"), paste(libnames,"(normalized and divided by gene length in kb)"), "Gene name", "Description")\n+ write.table(dataall, file=fileout, sep="\\t", quote=F, row.names=F)\n+\t\n+\t\t}\n+\t}\n+\telse{\n+\t\tprint("The files are not the same length")\n+\t}\n+}\n+\n+# Build a dataframe containing the files loaded into galaxy and their corresponding condition\n+\n+args <- commandArgs(trailingOnly = TRUE)\n+\n+Files=c()\n+Sample_name =c()\n+nbcells = (length(args)-3)\n+for (i in seq(1,nbcells,2)){\n+\tFiles = c(Files, args[3+i])\n+\tSample_name = c(Sample_name, args[4+i])\n+}\n+nfiles=nbcells/2\n+corresp = data.frame(Files=Files, Sample_name=Sample_name)\n+\n+# Take the informations given by the galaxy user to run the script\n+RNAseqDataAnnotation(args[1], args[2],args[3],corresp,nfiles)\n+\n+\n+\n+\n+\n+\n' |
| b |
| diff -r ea3f7a6504e0 -r f183f8648c5a RNAseqDataAnnotation/RNAseqDataAnnotation.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/RNAseqDataAnnotation/RNAseqDataAnnotation.xml Wed Dec 10 06:42:21 2014 -0500 |
| b |
| @@ -0,0 +1,57 @@ +<tool id="RNAseqDataAnnotation" name="RNAseqDataAnnotation" version="1.0.0"> + <description>tool for RNAseq Data Normalisation and Annotation</description> + <requirements> + </requirements> + + + <command interpreter="R --vanilla --silent --slave -f"> +./RNAseqDataAnnotation.R --args + $Species + $ensversion + $fileout + #for $i in $datafiles + ${i.file} + ${i.condition} + #end for + </command> + + <inputs> + <repeat name="datafiles" title="File" help="Upload the files to concatenate and annotate"> + <param name="file" label="Select file" type="data"/> + <param name="condition" label="Condition of the file" type="text"/> + </repeat> + <param name="Species" type="select" label="Select the specie for your data" help="If your specie of interest is not listed, your data will be normalized but no annotation will be added. Contact us if you want us to add your specie." > + <option value="hsapiens_gene_ensembl">Homo sapiens</option> + <option value="mmusculus_gene_ensembl">Mus musculus</option> + <option value="None">Other specie</option> + </param> + + <param name="ensversion" type="select" label="Select the version of Ensembl to use" > + <option value="oct2014">Version 77</option> + <option value="aug2014">Version 76</option> + <option value="feb2014">Version 75</option> + <option value="dec2013">Version 74</option> + <option value="sep2013">Version 73</option> + <option value="jun2013">Version 72</option> + <option value="apr2013">Version 71</option> + <option value="jan2013">Version 70</option> + <option value="oct2012">Version 69</option> + <option value="jul2012">Version 68</option> + <option value="may2012">Version 67</option> + </param> + </inputs> + + <outputs> + <data name="fileout" label="Data_normalization_annotation.txt" format="tabular"/> + </outputs> + <help> +**What it does** + This function concatenates htseq-count result files and normalizes data. + If the specie is known, data will be annotated using the specified Ensembl annotations. + +**How** + Normalization is performed using the method described in Genome Biology 2010;11(10):R106. + and implemented in the DESeq2 Bioconductor package. + Gene length correspond to the median of the size of all transcripts corresponding to this gene. + </help> + </tool> |