This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Function for barplot

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

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.

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.

2 answers

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.

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.

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.

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

Log in to answer this question.