Loop in R script

Question

[Seq225]

Hello,

I have this R script:

#!/usr/bin/env Rscript

library(ggplot2)

Pai <- read.csv("abc.txt", sep = "\t")

colnames(Pai) <- c("nt","a","c","g","t")

Pai[2:5] <- scale(Pai[2:5])

write.table(Pai, file = "PaiZscore.txt", row.names = FALSE)

pdf("abc.pdf")

ggplot() + 
    geom_line(aes(Pai$nt,Pai$t), color = 'red') + 

theme_classic()

dev.off()

It runs one file abc.txt and create one pdf abc.pdf. How can I modify this script so that it can run as a loop, use 200 files as input (abs1.txt, abc2.txt.........abc200.txt) and creates 200 pdfs (abs1.pdf, abc2.pdf.........abc200.pdf)?

The abc.txt files is generated using multiple custom pipelines, which incorporated bowtie, samtools, bedtools, and piPipes. I can use the aforementioned tools, but I am no expert in R. My goal is to make the heatmap of the final result in R.

Thanks a ton!!

[JC]

Check For and While loops in R https://cran.r-project.org/doc/contrib/Paradis-rdebuts_en.pdf

[cpad0112]

1. write a function for ggplot
2. Loop the function over the text files

[Ram]

Or,

1. Write a function that accepts a file name and outputs a ggplot
2. lapply that function over a list of file names. Avoid loops.

Answer 1

[bioinformatics2020]

Put this in the directory where you have the 200 files. Make sure you don't have any other text files because it uses the list.files() function to look for .txt files.

if(!require("tidyverse")) install.packages("tidyverse")
if(!require("magrittr")) install.packages("magrittr")
library(tidyverse) 
library(magrittr)
files <- list.files(pattern = "*\\.txt$")
Pai <- map(files, read_tsv)
Pai <- map(Pai, ~ bind_cols(.x[1], as_tibble(scale(.x[2:5]))))
iwalk(Pai, ~ write_tsv(.x, path = paste0("zscore_",.y,files[[.y]])))

plot.me <- function(x, x.cord, y.cord, filename, device = "pdf", height = 10, width = 10) { 
  ggplot(data = x, mapping = aes_string(x = x.cord, y = y.cord)) + 
    geom_point(color = "red") +
    theme_classic()
  ggsave(filename = paste(filename, device, sep = "."), device = device, height = height, width = width)
}
iwalk(Pai, ~plot.me(.x, x.cord = "nt", y.cord = "t", device = "pdf", filename = gsub(".txt","",files)[.y]))

[Ram]

Interesting. Why use map instead of lapply?

[bioinformatics2020]

The main reason here was because I wanted to use imap, which allows you to iterate over the list of objects and their indices, and so I just wanted to be consistent throughout the calls.

[Ram]

Thanks for the explanation - I've never heard of imap before - learned something new today!

[rpolicastro]

If you want to use imap to save a file, or any other purpose where you don't want to return a list, there's a handy alternative function called iwalk. There is also walk, pwalk, and walk2 that are alternatives to map, pmap, and map2.

[Seq225]

Thanks a lot!!

I am no expert in R. Even understanding a minor issue takes time for people like me.

Would you please help me fix this or direct me to somewhere where I can get some helpful insights? I have tried, but it seems to me that x needs to be defined as a function??

Error in .x[1] : object of type 'closure' is not subsettable
Calls: map -> .f -> bind_cols -> list2
Execution halted

Thanks!

[bioinformatics2020]

Can you send me what one of these files look like inside?

[Seq225]

Thanks for your time, and appreciate your help. I don't quite see the option of sending you the file. I am pasting it below.

Thanks!!

"nt" "a" "c" "g" "t"
2 1.25468332957509 0.158563316582844 -1.64006710598786 0.926245535065208
3 0.219803815359422 -0.514600017160217 -1.56107070301328 2.21008499426014
4 -0.610531627769753 -0.333275269517582 -0.450825748677823 1.3783121894941
5 -0.562435399693933 -1.41509014868781 1.52408432568676 -0.593822373913796
6 -0.936667286688371 0.143229299868879 1.08652268935237 -0.677621308205621
7 -0.390315752029343 0.86009458124674 0.534481632158826 -0.910368496624458
8 -0.341949320425008 0.678386483186255 -0.0672356500731077 -0.0358259093712009
9 0.472984521578764 0.0979939505626822 -0.443542392375202 0.0917374535121999
10 1.15281659932462 -0.368543507959701 0.189922853227132 -1.07473376536896
11 2.53625866532571 -0.486615436657231 -1.05366354726401 -0.615455926683496
12 1.06067719610071 -0.278456159765157 -0.794077258529559 0.261821937356734
13 0.215210355374653 -0.700524969817042 -0.448211210517908 0.853139046395189
14 0.836408267432568 -0.226703853355526 -0.50610455548746 0.0512056247597743
15 0.447585389898275 -0.11974908677562 -0.850289828967738 0.797936187603542
16 0.527025227281932 0.235233400152669 -0.563064136828472 0.112127698651457
17 1.31385790232006 0.581398827470427 -0.987926587814707 -0.270811051524834
18 -0.0979555341752051 -0.136233154743132 -0.784552869518439 1.22314739721488
19 1.66323106469104 -0.89411693083085 0.0464967598832076 -1.01256838384683
20 0.106588536911285 -1.38633886734912 1.10146290740903 -0.665436893427285
21 0.546479881335072 1.30976462138376 -0.92947298466803 -0.114900274667344
22 -1.05906948510605 0.983533415794158 0.156867620776775 0.127793374795032
23 -1.34251298651916 -1.34608707347497 2.4569141905994 -1.16275994560429
24 0.269251061077821 -0.00321055974948641 1.26076298243815 -1.92441020001336
25 0.984209597530732 -1.67998528742155 -0.24054217953035 0.504266925293022
26 -1.21740875281633 -0.451347198215112 0.209718642152205 1.13387790934911
27 -1.15877458712839 3.22996686439002 -1.76537818493809 1.32186602307201
28 0.690228158505499 0.590599237498806 0.673612412811462 -1.91520972354809
29 0.0282295136416919 -0.948552690165425 1.68375790616731 -1.65262315199863
30 -1.45194541556808 -0.681357448924586 1.6882399715843 -0.469740272395635
31 -1.50463510362867 1.45237097682364 -0.253988375781343 0.780778542303435
32 0.900176300162306 1.43358680634903 -1.10464704138235 -0.287471373772764
33 -1.09419594381311 1.61874505817016 -0.997637729551536 1.285312778737
34 -1.10878693435296 -0.463614411586284 -0.493405370139301 1.97808379041957
35 0.855592717957192 -1.40895654200222 0.93413246517445 -1.11725488633009
36 0.228180124743413 -0.163067683992571 -0.186944147252086 0.144702358569051
37 -1.21092386813195 0.728222037506642 0.167512526142144 0.418976021844667
38 -1.46896823786457 0.0255407215891979 1.0826008821125 -0.106197121254247
39 0.385168374811116 0.768473831380799 -0.0661151337188583 -0.764901503862685
40 -0.965849267768082 0.273568441937581 0.944030359636986 -0.545582037852627
41 -0.171721097460029 -1.16284557374309 0.447641614704497 0.316278811570116

[bioinformatics2020]

I made some minor adjustments! Please try now. :)

[Seq225]

Thank you so much! It works now!!! Much appreciation for your help!!