7

I have file1 like this:

CHR                     SNP     TEST   A1   A2                 GENO   O(HET)   E(HET)            P 
   0         AFFX-SNP-000541      ALL    0    0                0/0/0      nan      nan            1
   0         AFFX-SNP-000541      AFF    0    0                0/0/0      nan      nan           NA
   0         AFFX-SNP-000541    UNAFF    0    0                0/0/0      nan      nan           NA
   0         AFFX-SNP-002255      ALL    0    0                0/0/0      nan      nan            1
   0         AFFX-SNP-002255      AFF    0    0                0/0/0      nan      nan           NA
   0         AFFX-SNP-002255    UNAFF    0    0                0/0/0      nan      nan           NA
   1                 rs12103      ALL    C    T           55/250/317   0.4019   0.4113       0.5596
   1                 rs12103      AFF    C    T                0/0/0      nan      nan           NA
   1                 rs12103    UNAFF    C    T                0/0/0      nan      nan           NA
   1         rs12103_1247494      ALL    C    T           55/250/321   0.3994   0.4097       0.5581
   1         rs12103_1247494      AFF    C    T                0/0/0      nan      nan           NA
   1         rs12103_1247494    UNAFF    C    T                0/0/0      nan      nan           NA

And file2 like:

CHR                     SNP   A1   A2          MAF  NCHROBS
   0         AFFX-SNP-000541    0    0           NA        0
   0         AFFX-SNP-002255    0    0           NA        0
   1                 rs12103    C    T       0.2894     1244
   1         rs12103_1247494    C    T       0.2875     1252

I would like to merge file2 with file1 based on the SNP name and the TEST==ALL and keep the CHR, SNP, P, and MAF in the output file3. Does anybody know the best way for that in terminal (Unix) shell?

A desired output would be:

  CHR                     SNP  MAF        P
   0         AFFX-SNP-000541   NA         1
   0         AFFX-SNP-002255   NA         1
   1                 rs12103   0.2894     0.5596
   1         rs12103_1247494   0.2875     0.5581
4
  • Please give an example of your desired output.
    – terdon
    Feb 6, 2014 at 18:06
  • How large are the files? Is it possible that at the bottom of file1 there's again a line with 'AFFX-SNP-000541' for the SNP field?
    – brm
    Feb 6, 2014 at 18:18
  • No. file2's 'AFFX-SNP-000541' is exactly matching with file1's 'AFFX-SNP-000541'+'ALL'. Feb 6, 2014 at 18:20
  • Even if it is possible to do this with text files: wouldn't it at some point make sense to import the data into a relational database? Feb 7, 2014 at 9:44

4 Answers 4

6

With help from this answer

awk 'FNR==NR && FNR>1 {a[$2] = $5; next}
     FNR > 1 && ($2 in a) && $3 == "ALL" {
         print $1 "    " $2 "    "  a[$2] "    "  $9
     }' file2 file1

To get the header as well, just add this to the beginning of the script:

 BEGIN{print "CHR SNP MAF P"}

Explanation:

First of all, when two files are passed to awk, they are processed one after another. There are two variables important here: NR is the line number from the beginning of the awk command, and FNR is the line number from the beginning of the current file. That is, when the first file is processed (here file2), NR and FNR have the same value, which is the value of the line currently processed. But when awk pass to second file, FNR is reset to 1, so that NR and FNR are no longer the same. So that the test FNR==NR is a trick for knowing if the file processed is the first or not.

So let's see the code. The condition FNR==NR && FNR>1 tests whether we are processing the first file and not the first line. If it's the case, we store the value of fifth column (MAF) in an array indexed by the second one (SNP) and then the next statement says to pass to the following line.

When awk processes the second file (which is file1), the first test is false, so that awk tries the second test: FNR > 1 && ($2 in a) && $3 == "ALL", that is: not the first line of the file + second column value (SNP) exists in table a + third column value (TEST) is "ALL". If it is the case, then it prints column 1 (CHR) and two (SNP), gets the MAF value from the array with a[$2], and then prints column nine (P).

Adding a BEGIN{...} statement at the beginning adds a command which is run only before the first line is processed.

2
  • Well, if headers are known, one can add BEGIN{print "CHR SNP MAF P"} at the beginning.
    – JPG
    Feb 6, 2014 at 19:08
  • May you be able to give an explanation? Thanks for your great job! Feb 7, 2014 at 14:00
5

terdon mentions doing it with coreutils only, a humble suggestion:

(echo CHR SNP MAF P
 paste <(tail -n +2 file2) <(grep ALL file1) | 
   while read -r chr snp _ _ maf _ _ _ _ _ _ _ _ _ p; do 
     echo $chr $snp $maf $p; 
   done) | column -t

output:

CHR  SNP              MAF     P
0    AFFX-SNP-000541  NA      1
0    AFFX-SNP-002255  NA      1
1    rs12103          0.2894  0.5596
1    rs12103_1247494  0.2875  0.5581
4
  • +1 b/c paste is secretly my favorite Unix tool. It can even provide you world peace if you let it 8-)
    – slm
    Feb 6, 2014 at 22:02
  • @slm - "b/c", what's that?
    – grebneke
    Feb 7, 2014 at 8:00
  • @grebneke - b/c = because
    – slm
    Feb 7, 2014 at 8:01
  • 1
    @slm - thanks, that's sorted then. Only world peace to go and we're done.
    – grebneke
    Feb 7, 2014 at 8:03
4

Here's one way of doing it (copy/paste this directly into your terminal):

(
 awk -v OFS="\t" 'NR==1{print $1,$2,$5,"P"}' file2;
 awk '$3=="ALL"{print $2,$NF}' file1 | 
  while read snp p; do 
   awk -v snp="$snp" -v p=$p -v OFS="\t" '$2==snp{print $1,$2,$5,p}' file2;
  done
) > file3

The output looks like this:

CHR SNP MAF P
0   AFFX-SNP-000541 0   1
0   AFFX-SNP-002255 0   1
1   rs12103 C   0.5596
1   rs12103_1247494 C   0.5581

Notes:

  • I am 99% sure you can do this somehow using coreutils, perhaps a combination of sort,paste and join but I couldn't figure it out.
  • This is not efficient and can take a while if you are dealing with huge files.
  • This changes the field separators to tabs, remove all cases of -v OFS="\t" if you don't want that.
  • As you can see, the fields are no longer aligned. Since they are all separated by tabs, this is not a problem for any programs that will read the file, it is only an issue if you try to read it.

Explanation

  • The 1st awk line just prints the 1st, 2nd and 5th field of the header of file 2 followed by P. This is the header of the new file.
  • The second awk prints the 2nd and last ($NF) fields of each line of file1 where the 3rd field is ALL.
  • This is then passed through a shell while loop which reads the two fields into the variables $snp and $p.
  • These are then given to the last awk using the -v option and then, for each line of file2, if the second field is $snp, it prints fields 1,2,5 and the current value of $p.

Here's a perl script that does the same:

perl -le 'print "CHR\tSNP\tMAF\tP"; open($f1, "$ARGV[0]"); 
          while(<$f1>){
           s/^\s+//; @f=split(/\s+/); $k{$f[1]}=$f[$#f] if $f[2] eq "ALL"
          } 
          open($f2,"$ARGV[1]"); 
          while(<$f2>){
            s/^\s+//; @f=split(/\s+/); 
            next unless defined($k{$f[1]}); 
            print "$f[0]\t$f[1]\t$f[4]\t$k{$f[1]}"
          }' file1 file2
1

Here's another way to do it:

{ printf %s\\n "CHR SNP MAF P"; grep -F ALL file1 | sort -k2,2 | \
join -j2 -o 2.1,2.2,2.5,1.9 - <(sort -k2,2 file2); } | column -t

output:

CHR  SNP              MAF     P
0    AFFX-SNP-000541  NA      1
0    AFFX-SNP-002255  NA      1
1    rs12103          0.2894  0.5596
1    rs12103_1247494  0.2875  0.5581

How it works:

grep -F ALL file1 | sort -k2 prints the lines matching the string ALL, these are then sorted on 2nd field:

   0         AFFX-SNP-000541      ALL    0    0                0/0/0      nan      nan            1
   0         AFFX-SNP-002255      ALL    0    0                0/0/0      nan      nan            1
   1         rs12103_1247494      ALL    C    T           55/250/321   0.3994   0.4097       0.5581
   1                 rs12103      ALL    C    T           55/250/317   0.4019   0.4113       0.5596

this is then joined on 2nd field (-j2) with file2 (also sorted on 2nd field) using -o 2.1,2.2,2.5,1.9 to output 1st, 2nd and 5th field from file2 and 9th field from processed file1:

0 AFFX-SNP-000541 NA 1
0 AFFX-SNP-002255 NA 1
1 rs12103 0.2894 0.5596
1 rs12103_1247494 0.2875 0.5581

Since this was grouped {...} with printf %s\\n "CHR SNP MAF P" which prints the header, the entire output is then prettified with column -t.
Note that with this solution, lines order (from file2) is not preserved in the output. It just happens that your file2 was already sorted on 2nd field but if it wasn't e.g.:

CHR                     SNP   A1   A2          MAF  NCHROBS
   0         AFFX-SNP-000541    0    0           NA        0
   1                 rs12103    C    T       0.2894     1244
   0         AFFX-SNP-002255    0    0           NA        0
   1         rs12103_1247494    C    T       0.2875     1252

and if you wanted to preserve the lines order you could number the lines in file2 with nl -ba -nrz file2:

000001  CHR                     SNP   A1   A2          MAF  NCHROBS
000002     0         AFFX-SNP-000541    0    0           NA        0
000003     1                 rs12103    C    T       0.2894     1244
000004     0         AFFX-SNP-002255    0    0           NA        0
000005     1         rs12103_1247494    C    T       0.2875     1252

before sorting & joining and adjust the command to: join on 2nd field of (processed) file1 and 3rd field of (processed) file2 and output the line number field + the same fields as the first solution join -1 2 -2 3 -o 2.1,2.2,2.3,2.6,1.9 then re-sort the output with sort -k1n and delete the first field with cut -d' ' -f2- before piping the entire output to column -t:

{ printf %s\\n "CHR SNP MAF P"; grep -F ALL file1 | sort -k2,2 | join \
-1 2 -2 3 -o 2.1,2.2,2.3,2.6,1.9 - <(sort -k3,3 <(nl -ba -nrz file2)) | \
sort -k1n | cut -d' ' -f2-; } | column -t

output:

CHR  SNP              MAF     P
0    AFFX-SNP-000541  NA      1
1    rs12103          0.2894  0.5596
0    AFFX-SNP-002255  NA      1
1    rs12103_1247494  0.2875  0.5581

Note that both solutions assume there's only one occurrence of ALL on each line, namely the 3rd field (hence the grep -F ALL). If that's not the case you'll have to use a regex to filter only the relevant lines.

You must log in to answer this question.

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