3

I have multiple files in format: chr position value.

I want to combine them in format "chr", "position", "samp1", "samp2", "samp3", "samp4",........

For example:

Samp1:

chr position value
1   3774318 1
1   3774319 1
1   3775200 2
1   3775201 7
1   3775202 70
1   3775203 7
1   3775204 270
1   3775205 3
1   3775206 5

Samp 2:

chr position value
1   3775200 1
1   3775201 1
1   3775202 10
1   3775203 1
1   3775204 12
1   3775205 1
1   3775206 13
1   3775207 1
1   3775208 1
1   3775209 18

(and so on ...)

Desired output file: / I put random values in the output file

chr, position, value-samp1, value-samp2, value-samp3, value-samp4 
1 50204 2 17 5 2
1 50205 2 17 5 2
1 50206 2 18 5 2
1 50207 2 19 5 3
1 50208 3 19 5 3
1 50209 3 19 5 3

I tried join, merge, cat, but it does not work as I expected. I am a beginner. Do you have any ideas how it can be done?

(ed. note: here's an example of the operation described with actual data you provided.)

chr position    samp1   samp2
1   3774318     1       0
1   3774319     1       0
1   3775200     2       1
1   3775201     0       1
1   3775202     70      10
1   3775203     7       1
1   3775204     270     12
1   3775205     3       1
1   3775206     5       13
1   3775207     7       1
1   3775208     0       1
1   3775209     0       18
  • Please select your code blocks and format them with the {} button. – roaima Nov 8 at 16:05
  • Are all positions identical in all files or does their appearance vary? Are the files sorted by position? Does only position need to match or also chr? Sounds a bit, like you are trying to join a few files, doesn't it? – Fiximan Nov 8 at 16:17
  • Please can you provide a desired output that matches your sample inputs. – roaima Nov 8 at 16:18
  • 1
    @K7AAY well, join has the -a option and -o auto, too. I suggest OP reads man join at full attention and then adds file by file. For merging more files at once, I think one will have to resort to more elaborate languages rather than standard utils. It's just about my understanding of this sites philosophy of at least seeing some effort before helping where stuck rather than presenting full-blown solutions out of the box. – Fiximan Nov 8 at 17:07
  • 2
    is chr always 1? – rastafile Nov 8 at 17:47
2

Another awk approach. This one will print NA for cases where a particular variant isn't present in one of the input files:

awk '{ if(FNR==1){files[fnum++]=FILENAME}else{var[$1"\t"$2][FILENAME]=$3}} END{ for(v in var){for(file in files){if(! var[v][files[file]]){var[v][files[file]]="NA"}}}printf "chr\tposition "; for(i=1;i<=fnum;i++){printf "value-samp%s\t",i;} print "";for(v in var){ printf "%s ",v; for(file in var[v]){if(file in var[v]); else{var[v][file]="NA";}  printf "%s\t", var[v][file] } print ""}}' s1 s2

Or, if you're not into the whole brevity thing:

awk '{ 
        if(FNR==1){
            files[fnum++]=FILENAME
        }
        else{
            var[$1"\t"$2][FILENAME]=$3
        }
      } 
      END{ 
        for(v in var){
            for(file in files){
                if(! var[v][files[file]]){
                    var[v][files[file]]="NA"
                }
            }
        }
        printf "chr\tposition "; 
        for(i=1;i<=fnum;i++){
            printf "value-samp%s\t",i;
        } 
        print "";
        for(v in var){ 
            printf "%s ",v; 
            for(file in var[v]){
                if(file in var[v]); 
                else{
                    var[v][file]="NA";
                }
            printf "%s\t", var[v][file];
        } 
    print "";
    }
}' s1 s2

Change s1 and s2 to the actual file names you're using. Running the above on your example input returns:

chr position value-samp1    value-samp2 
1   3774318 1   NA  
1   3775200 2   1   
1   3774319 1   NA  
1   3775201 7   1   
1   3775202 70  10  
1   3775203 7   1   
1   3775204 270 12  
1   3775205 3   1   
1   3775206 5   13  
1   3775207 NA  1   
1   3775208 NA  1   
1   3775209 NA  18  
0

Taking this solution to a similar problem as foundation, we will need to fine tune it to get what you want. This is not a answer, just an outline, and uses some pseudocode (sic?) to outline what must be done.

A review of the Linux join command and its capabilities is essential to the successful solution. Please note sorting each input file into order on a specific field (here using the chromosome number, field 2) before joining is essential.

Because join can only join two files, we will need to step through the join multiple times, so will need some programming structure and control mechanism. We will also need to modify the data itself, as you specified blank (null) values should be replaced with 0.

1) Copy the first file of data to a different name, perhaps samp_0

2) Check each line of samp_0 for the value of field 3. If null, replace with 0 .

3) Initialize the value of x to match the number of the data file.

4) Create a for/next loop to use every data file consecutively. How many data files do you have? Well, to make this work without editing the script every times you run it, have the script run until it runs out of data files with a for/next loop.

5) Inside that for/next loop, do two things.

5A) Check each line of the next file for the value of field 3. If null, replace with 0 .

5B) perform a join including unpairable lines (with the -a option) on the chromosome number (field 2) of file samp_0 and file samp_x with automatic formatting (using the option -o auto) so it merges data lines for all chromosomes, not just those which have data in both file samp_0 and file samp_x . Write the output to file samp_0 .

6) If there's another data file, increment the value of x then go back and perform 3) again. If all data files have been joined, exit from the for/next loop, because you're finished.

0

I have ignored chr for the moment because if it is always 1 then it can be safely ignored, otherwise OP needs to explain how it fits in.

awk 'BEGIN {printf "position " }
        FNR>1{ 
        if (FNR==2) {nof+=1; printf ("%s%s " ,"Sam", nof )};
        pos[$2]=$2; data[$2, nof]=$3}
    END { printf "\n"; for (p in pos) {printf ("%s ", pos[p]); 
        for (d=1;d<=nof;d++) printf ("%s ", data[p,d]+0); print "\n"}
    }' file1 file2 | column -t

Output

position  Sam1  Sam2
3774318   1     0
3774319   1     0
3775200   2     1
3775201   7     1
3775202   70    10
3775203   7     1
3775204   270   12
3775205   3     1
3775206   5     13
3775207   0     1
3775208   0     1
3775209   0     18

Walkthrough

Print the first column header when you start

awk 'BEGIN {printf "position " }

Ignore all first lines as headers

FNR>1{ 

At the 2nd line of each file increase the count for the number of files and print a header

if (FNR==2) {nof+=1; printf ("%s%s " ,"Sam", nof )};

For every line after the header, put the position into an array pos; set another array data containing the sam value indexed by the pos and file (sample) number

pos[$2]=$2; data[$2, nof]=$3}

After reading all of the files in, print a newline on the titles then iterate over pos printing each position

END { printf "\n"; for (p in pos ) {printf ("%s ", pos[p]); 

Then iterate over the data array indexing by the p and sam/file number, adding 0 to any null data values so something prints even when there is no data; then print a newline

for (d=1;d<=nof;d++) printf ("%s ", data[p,d]+0); print "\n"}

Just the input files and piping the output through column to make it look nice

}' file1 file2 | column -t

Added chr on an assumption that it is another index

awk 'BEGIN {printf "chr position " }
        FNR>1{
        if (FNR==2) {nof+=1; printf ("%s%s " ,"Sam", nof )};
        chr[$1]=$1; pos[$2]=$2; data[$1, $2, nof]=$3}
    END { printf "\n"; for (c in chr) {for (p in pos ) {printf ("%s %s ", chr[c], pos[p]); 
        for (d=1;d<=nof;d++) printf ("%s ", data[c, p, d]+0); print "\n"}
    }}' file1 file2 | column -t
  • I think your assumption on chr is correct - there are many chromosomes, each with many more positions, at least in real life. Do you have any idea about performance/scaling? This must be some classical operation, this "de-normalisation". (normally you would have one long table with one extra field for the sample). – rastafile Nov 9 at 5:27
  • I deal with mass data a lot and TBH although scripting is a very useful hack for immediate problems with limited data, nothing beats a properly structured and indexed database. OP's question is a very, very routine SQL query in disguise. – bu5hman Nov 9 at 7:24
  • I doubt the chromosome will be constant. The OP should have clarified, but I work in the same field and so I would expect the chr values to range from 1 to 22 and also include Y, X and MT. – terdon Nov 9 at 13:52
  • @terdon With the associative arrays that should not be a problem – bu5hman Nov 9 at 19:25
0

I take another approach, more SQL style if I may say. I have only a basic idea and some parts. Idea is to use a "normalized" list i.e. first add the source-id to one large common table:

I have sample files ch.s1 to ch.sn; I take "chr" and "pos" as primary key.

. chr.sh >ch.tmp

for f in ch.s?
do 
  cat $f | while read l  
  do 
    echo "$l $f"
  done
done 

So now it is one big file with a 4th column, the source(file). For simplicity I used a tempfile.

From here, there is a "tee":

To get the list of sources for the comlumns you can do:

]# grep -o "ch.s.*" ch.tmp | sort | uniq -c
      3 ch.s1
      4 ch.s2
      2 ch.s3

(the -c could give valuable info; I use grep -o because I made a mess with delimiters (cut counts every space??)). This is just one way to get that list, but it is direct from the data. It is an investment.

This sorted list is the header, and can be used for filling the table.

And to prepare the reporting of that big table:

]# sort ch.tmp -k1,2 -k4 
1  16     A ch.s1
1  256    B ch.s1
1  333   dD ch.s2
1  4096   A ch.s1
1  4096  DD ch.s2
1  4096   A ch.s3
2  333   Dd ch.s2
2  777    x ch.s2
3 6666    F ch.s3

This sort is the only big operation - I count on it's efficiency in all situations.

After these two preparations, it is still a bit of looping, but only very "locally": for each data line you only loop through a fixed number of columns, putting the "value" or "null".

If you take this part:

1  4096   A ch.s1
1  4096  DD ch.s2
1  4096   A ch.s3

...it shows how the lines are ready to be printed out one by one, almost. No need for a full lookup (loop-in-a-loop).

With a elegant way to "tee" from ch.tmp (the unified table) to the "group by source" i.e. grep -o and the sort -k... this approach should be ready also for "bigger" data. Especially because the steps can be adjusted.

It simulates the reading into a SQL table, plus the reporting part as pivot table.


What would be a realistic scenario for benchmarking?

20 sample files, 50KB each, with a bit overlap as in the OP?

With perl, and also awk I think, you can do this approach very easy. With bash script I don't know, because of array performance; with bash you may use a tmp file, as in my outline.

0

Using Miller (https://github.com/johnkerl/miller) you can run

mlr --c2p --ifs ' '  --repifs \
sort -f chr,position,value \
then nest --implode --values --across-records -f value \
then nest --explode --across-fields --values -f value \
then unsparsify then clean-whitespace input0*.csv >output.csv

And obtain

chr position value_1 value_2
1   3774318  1       -
1   3774319  1       -
1   3775200  1       2
1   3775201  1       7
1   3775202  10      70
1   3775203  1       7
1   3775204  12      270
1   3775205  1       3
1   3775206  13      5
1   3775207  1       -
1   3775208  1       -
1   3775209  18      -

Note:

  • you must adapt the command to the names of your input files. In my example I have used input0*.csv
  • if you want a real CSV output change --c2p to --csv.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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