Take the 2-minute tour ×
Unix & Linux Stack Exchange is a question and answer site for users of Linux, FreeBSD and other Un*x-like operating systems.. It's 100% free, no registration required.

I need to sort a bed file randomly 10000 times and take the top 1000 rows each time. Currently, I am using the following code:

for i in {1..100}; do
    for j in {1..100}; do
        sort -R myfile.bed_sorted | tail -n 1000 > myfile.bed.$i.$j.bed
    done
done

It takes almost 6 hours to do this for each file. I have around 150 of them to be worked out. Is there a faster solution for this?

A sample of the data (myfile.bed_sorted) I have:

    chr1    111763899   111766405   peak1424    1000    .   3224.030    -1  -1
    chr1    144533459   144534584   peak1537    998 .   3219.260    -1  -1
    chr8    42149384    42151246    peak30658   998 .   3217.620    -1  -1
    chr2    70369299    70370655    peak16886   996 .   3211.600    -1  -1
    chr8    11348914    11352994    peak30334   990 .   3194.180    -1  -1
    chr21   26828820    26830352    peak19503   988 .   3187.820    -1  -1
    chr16   68789901    68791150    peak11894   988 .   3187.360    -1  -1
    chr6    11458964    11462245    peak26362   983 .   3169.750    -1  -1
    chr1    235113793   235117308   peak2894    982 .   3166.000    -1  -1
    chr6    16419968    16422194    peak26522   979 .   3158.520    -1  -1
    chr6    315344  321339  peak26159   978 .   3156.320    -1  -1
    chr1    111756584   111759633   peak1421    964 .   3110.520    -1  -1
    chrX    12995098    12997685    peak33121   961 .   3100.000    -1  -1
    chr9    37408601    37410262    peak32066   961 .   3100.000    -1  -1
    chr9    132648603   132651523   peak32810   961 .   3100.000    -1  -1
    chr8    146103178   146104943   peak31706   961 .   3100.000    -1  -1
    chr8    135611963   135614649   peak31592   961 .   3100.000    -1  -1
    chr8    128312253   128315935   peak31469   961 .   3100.000    -1  -1
    chr8    128221486   128223644   peak31465   961 .   3100.000    -1  -1
    chr8    101510621   101514237   peak31185   961 .   3100.000    -1  -1
    chr8    101504210   101508005   peak31184   961 .   3100.000    -1  -1
    chr7    8173062 8174642 peak28743   961 .   3100.000    -1  -1
    chr7    5563424 5570618 peak28669   961 .   3100.000    -1  -1
    chr7    55600455    55603724    peak29192   961 .   3100.000    -1  -1
    chr7    35767878    35770820    peak28976   961 .   3100.000    -1  -1
    chr7    28518260    28519837    peak28923   961 .   3100.000    -1  -1
    chr7    104652502   104654747   peak29684   961 .   3100.000    -1  -1
    chr6    6586316 6590136 peak26279   961 .   3100.000    -1  -1
    chr6    52362185    52364270    peak27366   961 .   3100.000    -1  -1
    chr6    407805  413348  peak26180   961 .   3100.000    -1  -1
    chr6    32936987    32941352    peak26978   961 .   3100.000    -1  -1
    chr6    226477  229964  peak26144   961 .   3100.000    -1  -1
    chr6    157017923   157020836   peak28371   961 .   3100.000    -1  -1
    chr6    137422769   137425128   peak28064   961 .   3100.000    -1  -1
    chr5    149789084   149793727   peak25705   961 .   3100.000    -1  -1
    chr5    149778033   149783125   peak25702   961 .   3100.000    -1  -1
    chr5    149183766   149185906   peak25695   961 .   3100.000    -1  -1
share|improve this question
1  
How big is your file and how strict is your notion of "random"? split can, err, split a file into pieces of 1000 lines each, so you'd get more files in a single call of sort. Also, have you checked if head is slightly faster than tail because it doesn't need to read through the entire file? –  Ulrich Schwarz Jun 30 at 11:38
    
@UlrichSchwarz: The sample file that I have pasted above contains around 33000 rows. In general, all my bed files will have more or less the same number of rows. Also for example: from a 33000 row file, I dont wish to get 33 subsets (1000 rows in each) in a single run. I only wish to take the top 1000 rows from each run. I will also be doing a tail of the same file. Just for sample, I used head here. –  biobudhan Jun 30 at 11:44
    
According to the man page sort -R uses a "random hash of keys". Creating the hash is a total waste of time and probably takes longer than anything else. It would be better to read the lines into an array and then shuffle that using indexes. Personally, I'd use perl for that; you could do it with bash but you'll need a function to generate random numbers. –  goldilocks Jun 30 at 12:19
    
@goldilocks: I am not a perl person! Could you please help me out? –  biobudhan Jun 30 at 12:37
6  
Try shuf instead of sort -R, it's considerably faster. Of course, doing it in memory (see Perl answer) will beat anything that requires re-reading the entire file in the shell. –  frostschutz Jun 30 at 13:30

3 Answers 3

up vote 14 down vote accepted

Assuming you have enough memory to slurp the file, you could try

perl -e 'use List::Util 'shuffle'; @k=shuffle(<>); print @k[0..999]' file.bed

Since you want to do this 10000 times, I would recommend integrating the repetition into the script and shuffling the indices instead of the array itself to speed things up:

$ time perl -e 'use List::Util 'shuffle'; 
            @l=<>; for $i (1..10000){
               open(my $fh, ">","file.$i.bed"); 
               @r=shuffle(0..$#l); 
               print $fh @l[@r[0..999]]
            }' file.bed

real    1m12.444s
user    1m8.536s
sys     0m3.244s

The above created 10000 files of 1000 lines each from a file that contained 37000 rows (your example file repeated 1000 times). As you can see, it took slightly more than three minutes on my system.

Explanation

  • use List::Util 'shuffle'; : this imports a Perl module that provides the shuffle() function which randomizes an array.
  • @l=<>; : load the input file (<>) into the array @l.
  • for $i (1..10000){} : run this 10000 times.
  • @r=shuffle(0..$#l); : $#l is the number of elements in @l so @r is now a randomized list of the index numbers of the array @l (the input file's lines).
  • open(my $fh, ">","file.$i.bed"); : open a file called file.$i.bed for writing. $i will take values from 1 to 10000.
  • print $fh @l[@r[0..999]] : take the first 1000 indices in the shuffled array and print the corresponding lines (elements of @l).

Another approach is to use shuf (thanks @frostschutz):

$ time for i in {1..10000}; do shuf -n 1000 file.bed > file.$i.abed; done

real    1m9.743s
user    0m23.732s
sys     0m31.764s
share|improve this answer
    
Wow!! That is awesome!! It worked in 2 mins :-) I have just one more question. How about also retrieving the last 1000 lines of the file? Because We need to know the length (number of lines) in the file to acheive this? Please help! –  biobudhan Jun 30 at 13:51
1  
@biobudhan do consider shuf as suggested by frostschutz: for i in {1..10000}; do shuf -n 1000 file.bed > file.$i.bed; done. That took ~1 minute on my system. As for the last 1000 lines, all you need is tail -n 1000. –  terdon Jun 30 at 13:56
1  
@biobudhan also see updated answer for a 3x faster perl version. –  terdon Jun 30 at 14:07
    
Yes, I tried it and it works faster now!! Thank you very much!!! :-) –  biobudhan Jun 30 at 14:16
    
Did you double check the output files of the perl version? It seems odd to me that it has so little sys time, which would be file I/O -- this should not be so totally different than the shuf one, which has ~30s sys. So I tested the perl one here (cut n' paste) and O_O it created 1000 files but all the files were empty... –  goldilocks Jun 30 at 15:08

If you want a benchmark to see how fast it can be done, copy paste this into 10kshuffle.cpp and compile g++ 10kshuffle.cpp -o 10kshuffle. You can then run it:

10kshuffle filename < inputfile

Where filename is a base path to use for the output files; they'll be named filename.0, filename.1, etc. and each contains the first 1000 lines of a shuffle. It writes the name of each file as it goes.

#include <cerrno>
#include <cstdlib>
#include <cstring>
#include <fcntl.h>
#include <fstream>
#include <iostream>
#include <string>
#include <sstream>
#include <unistd.h>
#include <vector>

using namespace std;

unsigned int randomSeed () {
    int in = open("/dev/urandom", O_RDONLY);
    if (!in) {
        cerr << strerror(errno);
        exit(1);
    }
    unsigned int x;
    read(in, &x, sizeof(x));
    close(in);
    return x;
}

int main (int argc, const char *argv[]) {
    char basepath[1024];
    strcpy(basepath,argv[1]);
    char *pathend = &basepath[strlen(basepath)];
// Read in.
    vector<char*> data;
    data.reserve(1<<16);
    while (!cin.eof()) {
        char *buf = new char[1024];
        cin.getline(buf,1023);
        data.push_back(buf);
    }

    srand(randomSeed());
    for (int n = 0; n < 10000; n++) {
        vector<char*> copy(data);
    // Fisher-Yates shuffle.
        int last = copy.size() - 1;
        for (int i = last; i > 0; i--) {
            int r = rand() % i;
            if (r == i) continue;
            char *t = copy[i];
            copy[i] = copy[r];
            copy[r] = t;
        }
    // Write out.
        sprintf(pathend, ".%d", n);
        ofstream file(basepath);
        for (int j = 0; j < 1000; j++) file << copy[j] << endl;
        cout << basepath << endl;
        file.close();
    }

    return 0;
}  

On a single 3.5 Ghz core, this runs in ~20 seconds:

   time ./10kshuffle tmp/test < data.txt
   tmp/test.0
   [...]
   tmp/test.9999
   real 19.95, user 9.46, sys 9.86, RSS 39408

data.txt was 37000 lines duplicated from the question. If you want the entire shuffle in the output file instead of the first 1000 lines, change line 54 to:

for (int j = 0; j < copy.size(); j++) file << copy[j] << endl; 
share|improve this answer
    
Thank you very much! Loved it :-) –  biobudhan Jul 1 at 7:12

So there is a Unix aspect to your question, but it is worth solving your fundamental problem first and then trying to find a Unix-y way to implement that solution.

You need to create 10,000 samples of size 1,000 each from a file with an unknown, large number of rows. It is possible to do this in one single pass of the file if you can hold 10,000 x 1,000 rows in memory. If you cannot hold that many rows in memory, you can still do it in one single pass if you know how many rows your file contains. If you do not know how many rows your file contains you need one additional pass to count the number of rows.

The algorithm, in the more difficult case when you don't know the number of rows, is to do the following for each sample (in parallel, keeping the samples in memory):

  • include the first 1,000 rows in the sample
  • for the n-th row (where n > 1000), include it with the probability 1000 / n and discard a random row from the rows you have already selected. (because of the likelihood of discarding some rows we need to hold the sample in memory until end of input)

An elegant way to implement the second step is to generate a random integer k in [1, n]. If k <= 1000 then include the row and replace the existing k-th row with it. Here is a more standard description of the algorithm: http://en.wikipedia.org/wiki/Reservoir_sampling

If you know the number of rows, R, then:

  • start with sample-size, s of 0
  • include n-th row with probability (1000 - s) / (R - n + 1) and output it immediately (and increment the sample size s)

How to do this on Unix? awk seems to be the answer per this post on the Internet (I cannot vouch for its correctness, but the code is there) https://news.ycombinator.com/item?id=4840043

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.