0

I have two codes to be done for all of my samples. 

cr1 = MEDIPS.seqCoverage(file = "1.bam", pattern = "CG", BSgenome = BSgenome, extend = extend, shift = shift, uniq = uniq)
MEDIPS.plotSeqCoverage(seqCoverageObj=cr1, type="pie", cov.level = c(0, 5, 10, 20, 30), main="cr1")

Then,

cr2 = MEDIPS.seqCoverage(file = "2.bam", pattern = "CG", BSgenome = BSgenome, extend = extend, shift = shift, uniq = uniq)    
MEDIPS.plotSeqCoverage(seqCoverageObj=cr2, type="pie", cov.level = c(0, 5, 10, 20, 30), main="cr2")

I don't want to repeat the code again and again for 100 times. I tried some for loops but they are not working because " File i.bam not found in ...". Well, I am not good at it at all. Does anyone can help me out?

So my code look like this:

for(i in 1:100){
  cr[i] = MEDIPS.seqCoverage(file = paste0(as.character(i),".bam"),
                             pattern = "CG", BSgenome = BSgenome, extend = extend, 
                             shift = shift, uniq = uniq)
  MEDIPS.plotSeqCoverage(seqCoverageObj=cr[i], type="pie", 
                         cov.level = c(0, 5, 10, 20, 30), main="cr",paste0(as.character(i)))
}

Reading bam alignment 1.bam

Total number of imported short reads: 17254741

Extending reads...

Creating GRange Object...

Keep at most one 1 read mapping to the same genomic location.

Number of remaining reads: 11148075

Loading chromosome lengths for BSgenome.Hsapiens.UCSC.hg19...

Get genomic sequence pattern positions...

Number of identified CG pattern: 26752702

Calculating sequence pattern coverage...

Error: object 'cr' not found

Jinyong Huang
  • 49
  • 5
  • 2
    What exactly did your loops look like? It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. – MrFlick Feb 19 '20 at 21:09
  • For sample #1 I need cr1 as a vector, "1.bam" as the bam file, "cr1" as the main title of the plot. For sample #2, I will need cr2, "2.bam", "cr2" etc – Jinyong Huang Feb 19 '20 at 21:15
  • 1
    @JinyongHuang ... It's so much easier if you show us the literal code you ran. Don't just describe it. Show us *exactly* what you tried. – Dason Feb 19 '20 at 21:26
  • Do you need all the `cr1`, `cr2`, objects after the loop? Or are you just using the for plotting and then don't need them any more. It would be easier to store them in a single list rather than 100 separate variables in your global environment if you plan to hang on to them. – MrFlick Feb 19 '20 at 21:58
  • @MrFlick The answer from ThomasIsCoding has solved my problem. Thank you anyway. – Jinyong Huang Feb 19 '20 at 22:03
  • OK, but it may have also just created a new problem if you need to keep using those values multiple times. Generally you want to gather related values in lists and avoid `get/assign` to make it easier to do things in an R-like way. – MrFlick Feb 19 '20 at 22:05
  • @MrFlick Each time I asked a question, I learned more than how to solve the problem. Thanks. – Jinyong Huang Feb 19 '20 at 22:12

2 Answers2

1

I am not sure if this is the thing you are looking after, but you can still have a try, where the number similar actions is assumed as 100, i.e., 

list2env(setNames(lapply(paste0(seq(100),".bam"),
                         function(v) MEDIPS.seqCoverage(file = v, pattern = "CG", BSgenome = BSgenome, extend = extend, shift = shift, uniq = uniq)),
                  paste0("cr",seq(100))),
         envir = .GlobalEnv)

sapply(paste0("cr",seq(100)),function(v) MEDIPS.plotSeqCoverage(seqCoverageObj=get(v), type="pie", cov.level = c(0, 5, 10, 20, 30), main=v))
ThomasIsCoding
  • 96,636
  • 9
  • 24
  • 81
0

Be sure that your working directory is set correctly. The " File i.bam not found in ..." means that file is not in the same directory you are currently in. 

getwd()

will show you what directory you are currently in. If the file is in downloads, you can use 

setwd()

to change your directory to downloads to access the file.

Update from comment

for(i in 1:100){
cr[i] = MEDIPS.seqCoverage(file = cat(paste0(as.character(i),".bam")),
pattern = "CG", BSgenome = BSgenome, extend = extend, 
shift = shift, uniq = uniq)

MEDIPS.plotSeqCoverage(seqCoverageObj=cr[i], type="pie", 
cov.level = c(0, 5, 10, 20, 30), main=cr[i])
}
htii
  • 21
  • 1
  • 6