0

I have around 500 GB free disk space on my SSD. I'm trying to run an operation on 10 gzipped files (each around 25GB in size). But I keep running out of storage when I do it in a parallel for loop since sort writes a lot of temporary files in the same directory and apparently doesn't clean up after itself.

I'm trying to randomly take a certain number of lines from these files.

bcftools view "${FILES[i]}".vcf.gz | awk '{printf("%f\t%s\n",rand(),$0);}' | sort -t $'\t' -T . -k1,1g | head -n "${SUBSET_COUNT[i]}" | cut -f 2- >> "${FILES[i]}"_"${SUBSET_COUNT[i]}"_subset.vcf &

This operation takes around 1 hour on each file (when I do it one by one), but I want to do it in parallel, since I need to repeat this on more batches of such files.

Improve this question
2
  • 1
    This seems like a very complex way of getting random lines. What do you really need here? Just N randomly selected lines from each file? Don't you also need to exclude the vcf headers (lines starting with #)? Why do you need to sort? What is ${SUBSET_COUNT[i]} and where is it set?
    – terdon ♦
    Jan 16 at 18:31
  • @terdon I wanted to come up with an efficient way to truly get a random subset of the variants, and I read somewhere on Stack Overflow that shuf is inefficient especially for large files, that is why I avoided shuf. Yes, I'm already excluding the headers in my original for loop, here just for the sake of simplicity I kept the command to a bare minimum.
    – paropunam
    Jan 16 at 20:32

1 Answer 1

Reset to default
2

I don't really understand why you're taking this approach and I suspect it wouldn't work anyway since rand() isn't getting reinitialized and is likely always giving the exact same output. Try it:

 $ for i in {1..10}; do awk -v i=$i 'BEGIN{print "Try "i",rand="rand()}'; done
Try 1,rand=0.924046
Try 2,rand=0.924046
Try 3,rand=0.924046
Try 4,rand=0.924046
Try 5,rand=0.924046
Try 6,rand=0.924046
Try 7,rand=0.924046
Try 8,rand=0.924046
Try 9,rand=0.924046
Try 10,rand=0.924046

In any case, the whole thing seems way too complicated, you can just use shuf instead:

bcftools view file.vcf.gz | shuf -n 100 > newfile.vcf

That will pick 100 random lines from its input. Of course, that won't be a valid VCF file, that requires headers. So if you want to generate valid VCFs, add the header to each file and then pick N random non-header lines:

zgrep '^#' file.vcf.gz > newfile.vcf &&
zgrep -v '^#' file.vcf.gz | shuf -n 100 >> newfile.vcf

Or, better, avoid writing uncompressed text:

zgrep '^#' file.vcf.gz | bgzip > newfile.vcf.gz &&
zgrep -v '^#' file.vcf.gz | shuf -n 100 | bgzip >> newfile.vcf.gz
Improve this answer
6
  • You're right! Even though rand is giving me different numbers on each line, I've now just noticed that every time I run it on the same file, it gives me the exact same set of numbers! Isn't shuf inefficient for a 30GB+ file? I read here: stackoverflow.com/a/15065490
    – paropunam
    Jan 16 at 20:36
  • 1
    @paropunam is it? You have linked to some random comment, that doesn't seem conclusive. In any case, your approach is several orders of magnitude more inefficient since you are sorting the files, and processing all lines multiple times. The answer linked to in the comment you mention is completely different to what you were doing here and does all of the work in awk. It also only samples from the first N lines of the file and it requires keeping the entire N lines in memory which is efficient in terms of speed, but not in terms of resources.
    – terdon ♦
    Jan 16 at 22:10
  • 1
    By the way, @paropunam, you might be interested in our sister site: Bioinformatics.
    – terdon ♦
    Jan 16 at 22:24
  • This is all very excellent info, I'll spend some more time understanding these commands. So far shuf does indeed work a lot faster than my terribly complex solution. Thank you for your help. And fantastic, I had no idea there's a bioinformatics SE. Thank you for pointing that out.
    – paropunam
    Jan 16 at 22:26
  • 1
    @paropunam yes, it's a younger site and not too active, but the community there knows much more about the various bioinformatics formats and existing tools, so when you ask there, you don't need to explain what a VCF is etc. And we need more users, so come on by! :)
    – terdon ♦
    Jan 16 at 22:29

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .