Function for barplot
2
0
Entering edit mode
10.2 years ago
dktarathym ▴ 40

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 • 3.2k views
ADD COMMENT
1
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY
1
Entering edit mode
10.2 years ago
Philippe ★ 1.9k

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 COMMENT
1
Entering edit mode
10.2 years ago
linus ▴ 360

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 COMMENT
2
Entering edit mode

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 REPLY
0
Entering edit mode

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 REPLY

Login before adding your answer.

Traffic: 1598 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6