Question: Function for barplot
0
gravatar for dktarathym
4.9 years ago by
dktarathym40
United States
dktarathym40 wrote:

Hi,

I wrote some codes. Explanations are below: 

# Loading required packages

library("VariantAnnotation")
library("TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")

 

# args will take different values (VCF files) while running the script from command line 

args <- commandArgs(trailingOnly = TRUE)

file1<- table(seqnames(readVcf(args[1],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
file2<- table(seqnames(readVcf(args[2],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
file3<- table(seqnames(readVcf(args[3],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
file4<- table(seqnames(readVcf(args[4],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
file5<- table(seqnames(readVcf(args[5],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
file6<- table(seqnames(readVcf(args[6],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))

all<- rbind(file1, file2, file3, file4, file5, file6)
all<- as.matrix(all)

pdf("Barplot_of_mutaton.pdf")

barplot(all, beside=TRUE, col= rainbow(nrow(all)),
        legend=rownames(all), cex.names=0.7, args.legend = list(x = "topright",bty = "n",
                                                                inset=c(-0.07, -0.14)))
dev.off()

 

This script was for 6 inputs. This worked properly.

Now, I want to make a function, that will take input, while I run the script through command line. Input could be 1 to n. But it should plot the barplot.

For example: if I input 2 files, than also it should plot a barplot and if I provide input of 10 files, than also script should work. I am not good at writing functions, but tried the following commands:

 

args <- commandArgs(trailingOnly = TRUE)
args<- NULL
for (i in args) {
  library("VariantAnnotation")
  library("TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")
  args[i]<- table(seqnames(readVcf(args[i],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
  all<- rbind(args[i])
  pdf("Barplot_of_mutaton.pdf")
  barplot(all, beside=TRUE, col= rainbow(nrow(all)),
          legend=rownames(all), cex.names=0.7, args.legend = list(x = "topright",bty = "n",
                                                                  inset=c(-0.07, -0.14)))
  dev.off()
}

 

 

Regards,

Deepak

R • 2.0k views
ADD COMMENTlink modified 4.4 years ago by Biostar ♦♦ 20 • written 4.9 years ago by dktarathym40
1

To make this clear, you want to make barplot for every file and don't know how to make such looping through files, yes?

PS.: I am asking this, because post title is about barplots, but it seems that your problem is looping.

ADD REPLYlink modified 4.9 years ago • written 4.9 years ago by PoGibas4.8k

Probably, sapply function is to be used. But, how to use it and where to use it.

 

How can I learn about making functions and loops?
​Suggestions are welcome.

ADD REPLYlink modified 4.9 years ago • written 4.9 years ago by dktarathym40
1
gravatar for Philippe
4.9 years ago by
Philippe1.9k
Barcelona, Spain.
Philippe1.9k wrote:

I see several problems there (in addition to not explicitly telling us what goes wrong or providing an error message):

  • you reassign args[i] within the loop, that might create some confusion. You better use a local variable to store the content of your file
  • Using for (i in args) will assign the content of args to i, no the numerical index. args[i] will likely generate an error. You should use for (i in 1:length(args)) insteas
  • You assign the arguments to args and then assign it NULL. your args variable will then be empty...
  • There is no need to reload each library within the loop. Do it outside the loop, it will save some computation time.
  • You can indeed use sapply to avoid the for loop. But you will still need to write a function. Both will work anyway.
ADD COMMENTlink written 4.9 years ago by Philippe1.9k
1
gravatar for linus
4.9 years ago by
linus330
Germany
linus330 wrote:

The idea of parsing the arguments looks good.

Your problem is that, each iteration of the for loop writes the file "Barplot_of_mutation.pdf"

This results in only one file containing the last barplot. If you want to have all barplots in one file put the pdf(...) line outside of the for loop and the dev.off() at the end of your code.

if you want to have a single file for each barplot try this one:

pdf(paste0("Barplot_of_mutaton_", i, ".pdf"))

 

I did not check each of your codelines, but put the library(...) statements outside of the for loop. Otherwise you will load them multiple times, which is unnecessary.

 

ADD COMMENTlink written 4.9 years ago by linus330
2

OP, can start

 pdf("Barplot_of_mutation.pdf")

Before loop and dev.off() after loop. Like this he will get all plots in one file on separate pages.

ADD REPLYlink written 4.9 years ago by PoGibas4.8k

My idea is that, I want to make a script, that would work to plot the mutation for vcf file(s). For example: if I have 1 vcf file, and I want to see the no. of mutations in it (for each chromosome), this script should print the barplot. In case, I have multiple files, say 5, it should plot mutations for each chromosome in a single barplot (for bar plot, command is:
barplot(all, beside=TRUE, col= rainbow(nrow(all)),
        legend=rownames(all), cex.names=0.7, args.legend = list(x = "topright",bty = "n",
                                                                inset=c(-0.07, -0.14))).

)

 

If I understood properly, you want me to do the following:
 

library("VariantAnnotation")
library("TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")

args <- commandArgs(trailingOnly = TRUE)

pdf("Barplot_of_mutaton.pdf")
for (i in 1:length(args)) {
  args[i]<- table(seqnames(readVcf(args[i],"TxDb.Scerevisiae.UCSC.sacCer3.sgdGene")))
  all<- rbind(args[i])
  barplot(all, beside=TRUE, col= rainbow(nrow(all)),
          legend=rownames(all), cex.names=0.7, args.legend = list(x = "topright",bty = "n",
                                                                  inset=c(-0.07, -0.14)))
}
dev.off()

 when I passed the files through command line, there was following error message:
 

Error in -0.01 * height : non-numeric argument to binary operator

Calls: barplot -> barplot.default

In addition: Warning message:

In args[i] <- table(seqnames(readVcf(args[i], "TxDb.Scerevisiae.UCSC.sacCer3.sgdGene"))) :

  number of items to replace is not a multiple of replacement length

Execution halted

ADD REPLYlink written 4.9 years ago by dktarathym40
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1944 users visited in the last hour