2

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..
3
  • 1
    It looks like that you don't need to form those myDuplicatedFiles* files explicitly. Dec 27, 2017 at 14:28
  • What do you mean by that? @WeijunZhou Dec 27, 2017 at 14:50
  • I mean that if the ranking is what you really need, you don't need intermediate files. By the way, AFAIK there is not a good way to typeset multiline contents in comment, let alone codes. Dec 27, 2017 at 21:27

4 Answers 4

3

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.

16
  • Thanks again @steeldriver for helping me out. it didn't quit work this time, I probably have a typo, so i used : tail -n +2 sampleUnique.txt | split -l1 --filter='{ head -n 1 standard.txt && cat <(tail -n +2 standard.txt) - | sort -g -r -k4,4 } > "$FILE" ' I received the flowing error ` bash: -c: line 1: sytax error: unexpected end of file split: with FILE=xaa, exit 1 from command: { head -n 1 standard.txt && cat <(tail -n +2 standard.txt) - sort -g -r -k4,4 } > "$FILE" ` Dec 27, 2017 at 16:06
  • i am not that good in formatting my comments yet sorry for that post, kinda new to stackexchange , basically it says i have a syntax error. the standard.txt file is myDuplicationFile and the sampleUnique.txt is --> myUniqueFile. i hope you can help me out with this. Dec 27, 2017 at 16:11
  • @OsmanAltun if you put the compound command on a single line, you will need to terminate it with a semicolon e.g. '{ command1 && command2; } > "$FILE"' Dec 27, 2017 at 16:17
  • it works just the way i needed it, 1 small question: so after I sort them for method1. I want to sort the sorted(method1) files again. I used something like: for f in xa*.txt; do sort -g -r -k,k $f>> $f.method2.txt; done , but if i look at the out put i expect a different result. I sort first on the phylop colum and than the GPS column so the highest phylop score = 10.263 and the highest GPS =5 should be placed number one but they dont, do you know maybe why? Is it something i do wrong in the code or? Sorry for bothering you again. Dec 28, 2017 at 9:43
  • 1
    @Wildcard you can find it gnu.org/software/coreutils/manual/html_node/… its there if you look good, but i don't think that my problem has something to do with that. Jan 4, 2018 at 11:43
3

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
1
  • 1
    Comments are not for extended discussion; this conversation has been moved to chat.
    – terdon
    Jan 6, 2018 at 13:12
2

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);
8
  • Note: sorting the arrays in perl will also take a lot of memory. If the dupefile is huge, you may want to pre-sort it into two different files, and read the @phylop and @GPS arrays from the two differently sorted files.
    – cas
    Jan 4, 2018 at 13:11
  • thank you for the comment, the reason I wanted to keep the files is to just look in them after wards so that the sorting went good and the rest of the information can quite be important, in this example I only use 5 columns , normally the files consist of 151 columns. So I want test your perl code for my real files, do I have to put all the column names , I think that's not so efficient , you might help me out with that, I'm completely new to perl! Jan 4, 2018 at 13:53
  • no, you don't need to list all the column names. As long as you know what column number(s) you're interested in. I only used named variables because there were only 5. you can split into an array instead - e.g. my @fields = split;...then (remembering that perl arrays start from 0, not 1) $fields[fieldnumber] contains the value you want. e.g. instead of $phylop, use $fields[3]
    – cas
    Jan 4, 2018 at 13:56
  • as for verifying the results, i tend to test code thoroughly against a small but representative sample until I'm certain that it's working perfectly. for large files, visual inspection is going to introduce huge risk of human error, tired eyes, etc anyway.
    – cas
    Jan 4, 2018 at 14:03
  • I didn't get the syntax correct for my @fields =split' and $fields[fieldnumber] ` , assume phylop is column 43 and GPS is 151 , can you maybe edit it in the code so I can see it. thank you very much! Jan 4, 2018 at 14:12
1

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.

1
  • thank you for the explanation as well for the code, i did go for the split version in the end. Dec 28, 2017 at 9:12

You must log in to answer this question.

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