Duplicate a file multiple times, write to duplicated files, sort the files, count position of specific lines after sorting

Question

Before I go to more in details, I want to point out that I already asked a certain part of this question -> You can find it here. I received some nice answers, but I need to do more so I'll repeat my question and add more details this time:

So I have a file with unique content that looks like this(lets call this myUniqueFile):

chromosoom  start    end       phylop   GPS
chr1    28745756    28745756    7.905   5   
chr1    31227215    31227215    10.263  5
chr1    47562402    47562402    2.322   4
chr1    64859630    64859630    1.714   3
chr1    70805699    70805699    1.913   2
chr1    89760653    89760653    -0.1    0
chr1    95630169    95630169    -1.651  -1

These are all different locations with different scores as you can see.

And I have another file that looks like this(lets call this one myDuplicationFile:

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5

So to start with I would like to add lines(except the header) from myUniqueFile to myDuplicationFile, but I want them to be added in way that myDublicationFile is duplicated for each new line that is added from myUniqueFile. So myDublicationFile keeps his standard content + 1 new line added from myUniqueFile. it should look something like this:

myDublicatedFile1.txt:

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    0.905   1    <- first line from `myUniquefile`


myDublicatedFile2.txt:

chromosoom  start    end       phylop   GPS
chr3    15540407    15540407    -1.391  -1
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    33093371    33093371    3.753   4
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    39597927    39597927    8.721   5
chr1    31227215    31227215    10.263  5    <- second line from `myUniquefile`

so for each new line added a new file is created like myDublicatedFile3,4,5 etc..

After I have these myDublicatedFiles with the new added content, I would like to sort these files for a specific columns from high to low , (for the phylop column) I do this with for f in myDublicatedFile* ; do sort -g -r -k 3 $f >> $f.method1.txt so this looks something like this:

myDublicatedFile1.method1.txt:

chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    7.905   5 <- count 2
chr3    33093371    33093371    3.753   4
chr3    30648039    30648039    2.214   3
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    31663820    31663820    0.713   3
chr3    15540407    15540407    -1.391  -1
chromosoom  start    end       phylop   GPS

So after I sort these files, I would like to know the position of the lines I added after the sorting. It seems logical to me to do something with "grep" and use a "count".

So for myDublicatedFile1.method1.txt is this count/rank 2 since the added line from myUniquefile ended up in the second place in the file.

After calculating the count/rank for the phlop(method1) column i would like to do a sort for the GPS(method2) column and than calculate the ranks of the added lines again. myDublicatedFile1.method1.method2.txt should look something like this:

chr3    39597927    39597927    8.721   5
chr1    28745756    28745756    7.905   5 
chr3    33093371    33093371    3.753   4
chr3    30648039    30648039    2.214   3
chr3    31663820    31663820    0.713   3
chr3    37050398    37050398    1.650   2
chr3    38053456    38053456    1.1     1
chr3    15540407    15540407    -1.391  -1
chromosoom  start    end       phylop   GPS

Its easy if the counts/rankings are written in different file so I can use these later for statistics. So the most important files are these count since I will end up use these :)

Something like:

countsForMethod1.txt:

29
3
5
6
50
etc.

countsForMethod2.txt:

7
3
21
45
etc..

Answer 1

Assuming you have the version of split from GNU coreutils, and a shell like bash, ksh or zsh is available (for the process substitution feature being used here), you could modify the previous accepted answer to handle the header lines and sorting e.g.

tail -n +2 myUniqueFile | SHELL=$(command -v bash) split -l1 --filter='{ 
  head -n 1 myDuplicationFile &&
    sort -g -r -k4,4 <(tail -n +2 myDuplicationFile) -
  } > "$FILE"'

You could then use a simple awk one-liner to find the positions of the myUniqueFile entries in the output files:

awk 'FNR==NR && NR>1 {a[$0]++; next} ($0 in a) {print FILENAME, FNR}' myUniqueFile xa?
xaa 3
xab 2
xac 4
xad 5
xae 5
xaf 8
xag 9

Rinse and repeat for the other methods / sort orders.

Answer 2

This script calculates ranks without creating temporary files (almost, one file will created - the sorted_file). Also, it sorts myDuplicationFile one time for the each method, then use it further.

#!/bin/bash

rank_determination() {
    # Sorts the "myDuplicationFile" one time
    # The "sorted_file" will be used further.
    ###
    tail -n +2 myDuplicationFile | sort -g -r -k "$1","$1" > sorted_file

    # gawk iterates through "myUniqueFile" line by line (except the first line).
    gawk -v field_number="$1" '
    NR != 1 {
        # Stores the needed value for the each line
        ###
        search_value=$field_number
        cnt=1

        # then, it checks the specified column in the "sorted_file"
        # line by line for the value, which is less than 
        # the "search_value" from the "myUniqueFile".
        ###
        while((getline < "sorted_file") > 0) {
            if($field_number < search_value)
                break
            cnt++
        }

        print cnt
        # closing is needed for reading the file from the beginning
        # each time. Else, "getline" will read line by line consistently.
        ###
        close("sorted_file")
    }' myUniqueFile
}

# I create a function, which takes
# the number argument, which means the column number:
# "4" for "phylop" column, "5" for the "GPS" column.
#
# The function creates output, which you can redirect
# to the needed file.
# Call this function multiple times with different arguments
# for the each needed column.
rank_determination 4 > method1.txt
rank_determination 5 > method2.txt

Output

tail -n +1 -- method*
==> method1.txt <==
2
1
3
4
4
7
8

==> method2.txt <==
2
2
3
5
6
7
8

Answer 3

I agree with what @WeijunZhou said in his comment, you don't need to make all those temporary files to do this.

The following perl script will calculate the counts for your Method 1 (phylops) and Method 2 (GPS) sorting, in one pass through both files.

It works by keeping a sorted list (array) of the phylop and GPS values from the duplicate file, and then (for each line of the unique file) calculating where the phylop and GPS value would have sorted into their respective sorted arrays.

#!/usr/bin/perl

use strict;

# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';

# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
my @phylop=();
my @GPS=();

open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
  chomp;
  next if (m/^chromosoom/);

  my($chr,$start,$end,$phylop,$GPS) = split;
  push @phylop, $phylop + 0; # add 0 to make sure we only ever store a number
  push @GPS, $GPS + 0;
};
close(DUPE);

# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$a <=> $b} @phylop;
@GPS = sort {$a <=> $b} @GPS;

print "Method1\tMethod2\n";

# Now find out where the phylop and GPS value from each line of uniqfile
# would have ended up if we had sorted it into dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
  next if (m/^chromosoom/);
  chomp;

  my $phylop_sort_line=1;
  my $GPS_sort_line=1;

  my($chr,$start,$end,$phylop,$GPS) = split;

  for my $i (0..@phylop-1) {
    $phylop_sort_line++ if ($phylop < $phylop[$i]);
    $GPS_sort_line++ if ($GPS < $GPS[$i]);
  };

  #printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
  printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;  
};
close(UNIQ);

when run against the sample data you provided above, the output is:

$ ./counts-for-methods.pl
Method1 Method2
2       1
1       1
3       2
4       3
4       5
7       7
8       7

The script completely ignores the header line in both files, so these line numbers may be off by one if your current algorithm counts them.

It also assumes that the values from the unique file will always sort immediately before an identical value in the duplicate file. If this is not what you want, change the < comparisons in the for my $i (0..@phylop) loop to <=.

If you need the values for Method 1 and Method 2 separately, you can easily extract them with awk. Or the perl script can easily be modified to open two output files, one for each method, and print the respective values to each file.

---

here's a version for dealing with 151 fields in the input lines. I don't have such an input file, so I tested it with the '5 fields version' that has been commented out in the code. The output was the same as the version above.

#!/usr/bin/perl

use strict;

# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';

my @phylop=();
my @GPS=();

# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
  chomp;
  next if (m/^chromosoom/);

  my @fields = split;

# 151 fields version:
  push @phylop, $fields[42]+0;
  push @GPS, $fields[150]+0;

# 5 fields version:
#  push @phylop, $fields[3]+0;
#  push @GPS, $fields[4]+0;

};
close(DUPE);

# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$b <=> $a} @phylop;
@GPS = sort {$b <=> $a} @GPS;

print "Method1\tMethod2\n";

# Now find out where the phylop and GPS from each line of uniqfile
# would have ended up if we had sorted it into the dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
  next if (m/^chromosoom/);
  chomp;

  my $phylop_sort_line=1;
  my $GPS_sort_line=1;

  my @fields = split;

  for my $i (0..@phylop-1) {

# 151 fields version:
    $phylop_sort_line++ if ($fields[42] < $phylop[$i]);
    $GPS_sort_line++ if ($fields[150] < $GPS[$i]);

# 5 fields version:
#    $phylop_sort_line++ if ($fields[3] < $phylop[$i]);
#    $GPS_sort_line++ if ($fields[4] < $GPS[$i]);
  };

  #printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
  printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;

};
close(UNIQ);

Answer 4

Just to save some typing, I will call the myUniqueFile sample and myDuplicationFile standard.

#!/bin/bash                                                                     

(
while read line; do
  echo $line|cat standard -|tail -n +2|sort -g -r -k 4|awk '/^chr1/{print FNR}' >> countsForMethod1.txt
  echo $line|cat standard -|tail -n +2|sort -g -r -k 5|awk '/^chr1/{print FNR}' >> countsForMethod2.txt
done
) <(tail -n +2 sample)

Explanation: The whole while loop is wrapped by a pair of parentheses to make bash treat it as a single command. This command takes the file sample as its input, with the header line stripped using the tail command. It is then consumed by the read command one line a time. That means inside the loop, $line is one line from the file sample. The variable is echoed and pipelined to cat to generate the myDuplicated* files, only that they are generated "on the fly" and never written to disk. The header lines are dropped by the tail before the files are sorted. awk is then used to find out which line the sample is on.

Edit: I think the split has its advantages while this answer eliminates the need for intermediate files.