5

I have two text files: 1st file is a Tab delimited file which looks like this:

chrom   pos ref alt a1  a2  a3  a4
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd

2nd file looks like this:

a1
a4

I want to extract those columns in the 1st file which are present in the 2nd file along with first 4 columns of the first file. So in the above case, the output will look like this:

chrom   pos ref alt a1  a4
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd

I want to do this in shell. How can I do this? I have a bigger file than shown here, so I have many columns in 1st file

cut -f 1-4,$(grep -Fwf file2.txt <(head -1 file1.txt)) file1.txt
5
  • What did you try so far? Mar 22 at 4:31
  • I added in the post. My idea is to get the index of those columns which are in my file2 and then use cut command on the first file to get the desired output. I am just unable to get the correct index numbers so that I can use it with cut Mar 22 at 4:37
  • Would join do this?
    – user10489
    Mar 22 at 7:22
  • @user10489 No, it would not. The join utility uses numerical indexing of fields and compares the values of those fields from two files. In this case, the field indexes in the first file are unknown as the field names are given by the names in the second file. If the second file had contained only the field indexes (e.g., 5 and 8), then the ordinary cut command could have been used as in cut -f "1-4,$(paste -s -d , - <file2)" file1
    – Kusalananda
    Mar 22 at 9:18
  • 2
    Are these vcf files? Have you looked into dedicated tools like bcftools?
    – terdon
    Mar 22 at 15:48

6 Answers 6

4

If perl is okay:

$ perl -F'\t' -lane 'if(!$#ARGV){ $h{$_}=1; close ARGV if eof; next }
                     @i = grep { exists $h{$F[$_]} } 4..$#F if $.==1;
                     print join "\t", @F[0..3, @i]' f2.txt f1.tsv
chrom   pos ref alt a1  a4
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd

The hash variable will use the lines from the second file as keys.

Then, grep is used to get the index from the header line of the TSV file by testing field names against the hash keys.

Finally, the first 4 columns and the filtered index values are used for printing.

1
  • 1
    You could leave the exists out, also. @i = grep { $h{$F[$_]} } 4..$#F if $.==1; works fine.
    – terdon
    Mar 22 at 15:59
4

Using the standard paste command, we can create a comma-delimited list of the field names in your second file (here called file2):

$ paste -s -d , - <file2
a1,a4

We can use this with Miller (mlr), a utility for processing structured documents (for example TSV files), and its cut subcommand to extract the first four fields along with the fields in file2 from the first file (here called file1):

$ mlr --tsv cut -f "chrom,pos,ref,alt,$(paste -s -d , - <file2)" file1
chrom   pos     ref     alt     a1      a4
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd

The cut subcommand of Miller allows you to use the field names to select what fields to extract, and here we select the four first fields by name along with the names from the file2 file.

Since the cut subcommand also accepts a list of field names that end with a trailing comma, you could opt for using tr '\n' ',' <file2 in place of the paste command mentioned earlier.


If we had been given the field indices that we would want to retain, e.g. as the file file3 like this:

5
8

... then we could have used the standard cut command like below to extract the wanted data:

$ cut -f "1-4,$(paste -s -d , - <file3)" file1
chrom   pos     ref     alt     a1      a4
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd

We can create the numerical field index list from the original files with a somewhat convoluted pipeline:

$ head -n 1 file1 | tr '\t' '\n' | grep -xF -f file2 -n | cut -d : -f 1 | paste -s -d , -
5,8

This extracts the header line from the data and replaces the tabs with newlines to get each field name on a separate line. The grep then outputs the line numbers from this that correspond to the field names in your second file. These are outputted on the form 5:a1, where the number before the : is the line number and the text at the end is the matching field name.

The number is isolated using cut, and the paste command is then used to get all the field indices into a comma-delimited list.

The full command, which emulates what the mlr command at the top of this answer does, would therefore look something like this:

cut -f "1-4,$(
    head -n 1 file1 | tr '\t' '\n' |
    grep -xF -f file2 -n | cut -d : -f 1 |
    paste -s -d , -
)" file1
1
  • 1
    See also the -o option of mlr cut so the columns are displayed in the order specified (would not make a difference with the OP's sample) Mar 22 at 16:18
3

Using any awk:

$ cat tst.awk
BEGIN { FS=OFS="\t" }
NR == FNR {
    a[$1]
    next
}
FNR == 1 {
    for ( inFldNr=1; inFldNr<=NF; inFldNr++ ) {
        if ( (inFldNr <= 4) || ($inFldNr in a) ) {
            out2in[++numOutFlds] = inFldNr
        }
    }
}
{
    for ( outFldNr=1; outFldNr<=numOutFlds; outFldNr++ ) {
        inFldNr = out2in[outFldNr]
        printf "%s%s", $inFldNr, (outFldNr<numOutFlds ? OFS : ORS)
    }
}

$ awk -f tst.awk file2 file1
chrom   pos     ref     alt     a1      a4
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd
10      12345   C       T       aa      dd

Note that, for efficiency, when reading each line of file1 the above only loops as many times as the number of output fields to be printed so if your input had 1000 fields but you only wanted to print 10, this would iterate 10 times per line rather than 1000.

0

I had a very similar problem recently and it is worth mentioning how I solved it as a complement to the excellent answers people have provided to this post.

  • move the static columns to another file (called fixed_cols) and transpose the remaining columns
  • grep -f filters #filters is the second file (see below)
  • transpose again and we have the desired (filtered) columns
  • paste fixed_cols and the filtered columns and we are done

My assumptions: I used this transpose function:

 transpose () 
{ 
    awk '{ for (i=1; i<=NF; i++) RtoC[i]= (i in RtoC?RtoC[i] OFS :"") $i; } 
    END{ for (i=1; i<=NF; i++) print RtoC[i] }' $1
}
# space separated columns only

$ cat filters
a1
a4

And the main commands:

$ awk '{print $1,$2,$3,$4 >"fixed_cols"; print $5,$6,$7,$8}' input_file | transpose | grep -f filters | transpose >filtered
$ paste  fixed_cols filtered

for a tab-separated result you could use tr, For example:

$ paste  fixed_cols filtered | tr ' ' '\t'
$ paste  output filtered | tr ' ' '\t'
chrom   pos ref alt a1  a4
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd

(don't forget to remove temporary files fixed_cols and filtered) Be careful when using tr so you don't convert spaces that may be inside each column. Using tr is just an example. there are many other ways to format your output.

0

Using Raku (formerly known as Perl_6)

~$ raku -e 'my ($header,@a) = lines.map: *.split(/ \s+ /); 
            $header .= list;  my @ind = <a1 a4>; 
            my @col = (0...3, $header.grep( / @ind /, :k ).Slip); 
            put $header[@col].join("\t"); 
            say $_.join("\t") for @a.map: *.[@col];'  data.csv

Above is an answer coded in Raku, a member of the Perl-family of programming languages. Like Perl, Raku is well-positioned to take on text-wrangling jobs in bioinformatics, etc.

First, the data is read in using Raku's lines routine. These lines are split on whitespace. Saving the data into $header scalar and @a array, Raku understands that $header takes the first line, and the next statement upgrades $header from a lightweight Seq to a list. For now we take the desired column names inline, and store them in @ind array.

We grep the header for the desired column names, returning :k their index positions (not values). Then a @col array is created with the desired 0...3 first four columns, along with the grep results, Slipped together (i.e. flattened) to create a simple index.

If we were working on rows we'd be using the @array[@index] idiom (which works for the single-line header). Because we're working on columns we need to map into the array elements, e.g. @a.map: *.[@col]. Finally the data is re-joined on tabs and output.

Sample Input:

chrom   pos ref alt a1  a2  a3  a4
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd
10  12345   C   T   aa  bb  cc  dd

Sample Output:

chrom   pos ref alt a1  a4
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd

Upgrading the code to work on two files is simple. Just use a path to create @ind index, which helps you conceptually separate column-index files from data files:

my @ind = "/path/to/index.csv".IO.lines;

If you really prefer taking data off the command line, Raku has the @*ARGS array for that. Use Raku's IO routine and make sure to input the files in your desired order:

~$ raku -e 'my ($header,@a) = @*ARGS[0].IO.lines.map: *.split(/ \s+ /);
            $header .= list;  my @ind = @*ARGS[1].IO.lines;
            my @col = (0...3, $header.grep( / @ind /, :k ).Slip);
            put $header[@col].join("\t"); 
            say $_.join("\t") for @a.map: *.[@col];'   data.csv  index.csv
chrom   pos ref alt a1  a4
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd
10  12345   C   T   aa  dd

https://docs.raku.org/routine/grep
https://raku.org

-1

One way this can work is like this:

cut -f 1-4,$(grep -wFf file2.txt <(awk -F'\t' 'NR==1 {for (i=1; i<=NF; i++) print i "\t" $i}' file1.txt) | cut -f 1 | tr "\n" "," | sed 's/,$//') file1.txt
3
  • 1
    If you add an exit after the print statement in your awk script then it won't have to read the rest of file1.txt discarding every line and so it'd be much more efficient.
    – Ed Morton
    Mar 22 at 8:40
  • tr "\n" "," changes it's input into something that is no longer a valid test file per POSIX (no terminating newline) so YMMV regarding what sed or any other text-processing tool does with it as that's then undefined behavior. You'll probably be OK with GNU sed, idk about others.
    – Ed Morton
    Mar 22 at 8:41
  • Using paste -s -d , - in place of the tr command would solve the issue that @EdMorton points out (and also makes your sed call at the end superfluous). I have not looked at the rest of your code though.
    – Kusalananda
    Mar 22 at 9:23

You must log in to answer this question.

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