4

I have a directory ballgown, in which there are around 1000 subdirectories as sample names. Each subdirectory has a file t_data.ctab. The filename is the same in all subdirectories.

ballgown
      |_______TCGA-A2-A0T3-01A
                   |___________ t_data.ctab
      |_______TCGA-A7-A4SA-01A
                   |___________ t_data.ctab
      |_______TCGA-A7-A6VW-01A
                   |___________ t_data.ctab

Like above ballgown has 1000 subdirectories. The t_data.ctab file in all those 1000 subdirectories looks like below with columns:

t_id    chr     strand  start   end     t_name  num_exons       length  gene_id gene_name       cov     FPKM
1       1       -       10060   10614   MSTRG.1.1       1       555     MSTRG.1 .       0.000000        0.000000
2       1       +       11140   30023   MSTRG.10.1      12      3981    MSTRG.10        .       2.052715        0.284182
3       1       -       11694   29342   MSTRG.11.1      8       6356    MSTRG.11        .       0.557588        0.077194
4       1       +       11869   14409   ENST00000456328.2       3       1657    MSTRG.10        DDX11L1 0.000000        0.000000
5       1       +       11937   29347   MSTRG.10.3      12      3544    MSTRG.10        .       0.000000        0.000000
6       1       -       11959   30203   MSTRG.11.2      11      4547    MSTRG.11        .       0.369929        0.051214
7       1       +       12010   13670   ENST00000450305.2       6       632     MSTRG.10        DDX11L1 0.000000        0.000000
8       1       +       12108   26994   MSTRG.10.5      10      5569    MSTRG.10        .       0.057091        0.007904
9       1       +       12804   199997  MSTRG.10.6      12      3567    MSTRG.10        .       0.000000        0.000000
10      1       +       13010   31097   MSTRG.10.7      12      4375    MSTRG.10        .       0.000000        0.000000
11      1       -       13068   26832   MSTRG.11.3      9       5457    MSTRG.11        .       0.995280        0.137788

From all the t_data.ctab files I want to extract only t_name and FPKM column and create a new file. In the new file the FPKM column should be the sample name. It should look like below:

t_name         TCGA-A2-A0T3-01A TCGA-A7-A4SA-01A    TCGA-A7-A6VW-01A
MSTRG.1.1              0            0.028181                 0
MSTRG.10.1         0.284182         0.002072             0.046302
MSTRG.11.1         0.077194         0.685535             0.105849
ENST00000456328.2      0            0.307315             0.038961
MSTRG.10.3             0            0.446015             0.009946
MSTRG.11.2         0.051214         0.053577             0.036081
ENST00000450305.2      0            0.110438             0.040319
MSTRG.10.5         0.007904             0                1.430825
MSTRG.10.6             0                0                0.221105
MSTRG.10.7             0            0.199354                 0
MSTRG.11.3         0.137788         0.004792                 0

If it is two or three files I can use cut -f6,12 on each file and then join them. But I have around 1000 files now.

  • Is the t_name column exactly the same in all files? – simlev Aug 29 '18 at 9:51
  • @simlev Yes, it is. – beginner Aug 29 '18 at 9:53
  • You could run cut -f6 on any file to create the first column, then find all files and use awk or perl on them, to add a subsequent column for each. – simlev Aug 29 '18 at 12:33
  • I'd recommend you output in CSV format: with 1000 columns, you'll want to use a spreadsheet program to view the data – glenn jackman Aug 29 '18 at 15:32
3

Try this simple way:

first do:

awk 'FNR==1 { print substr(FILENAME,1,16) >substr(FILENAME,1,16)".tmp" }
     FNR >1 { print $12 > substr(FILENAME,1,16)".tmp" }
     NR==FNR{ print $6  >"first_column.tmp" }' TCGA-A*/t_data.ctab

then paste them together with comma delimited file (remove -d, if you want to have Tab instead):

paste -d, *.tmp
t_name,TCGA-A2-A0T3-01A,TCGA-A7-A4SA-01A,TCGA-A7-A6VW-01A
MSTRG.1.1,0.000000,0.00000,0.0000
MSTRG.10.1,0.284182,0.28418,0.2841
MSTRG.11.1,0.077194,0.07719,0.0771
ENST00000456328.2,0.000000,0.00000,0.0000
MSTRG.10.3,0.000000,0.00000,0.0000
MSTRG.11.2,0.051214,0.05121,0.0512
ENST00000450305.2,0.000000,0.00000,0.0000
MSTRG.10.5,0.007904,0.00790,0.0079
MSTRG.10.6,0.000000,0.00000,0.0000
MSTRG.10.7,0.000000,0.00000,0.0000
MSTRG.11.3,0.137788,0.13778,0.1377
1

Would you be satisfied with csv output?

find ballgown -name t_data.ctab | awk ' {
  F=$0
  print F " started"
  split(F,P,"/")
  FN= P[2]
  TF[FN]=1
  getline < F
  while ((getline < F) > 0) {
    TN[$6]=1
    TV[FN ":" $6] = $NF
  }
  close(F)
  print f " done"
}
END {
  printf("tname")
  for (F in TF) {
    printf(", %s",F)
  }
  print ""
  for (N in TN) {
    printf("%s",N)
    for (F in TF) {
      printf(", %s",TV[F ":" N])
    }
    print ""
  }
}
'
  • I wrote this code in a shell script "eg.sh". Then I gave sh eg.sh > eg.csv; Is this right? – beginner Aug 29 '18 at 12:31
  • If you made the script executable and put it somewhere on $PATH, then yes. – Gerard H. Pille Aug 29 '18 at 12:32
  • In the starting of shell script I gave #!/bin/bash; and then I gave the code you posted. As I said before on the command line I gave sh eg.sh > eg.csv; looks like it is running. but its been one hour but didn't finish yet anything. – beginner Aug 29 '18 at 14:08
  • It depends on the size and number of those files. To see if it advances, I'd add a "print F" after the "F=$0" and "close(F)" statements, and restart the script. – Gerard H. Pille Aug 29 '18 at 14:27
  • could you please make this changes in the above code. thank you – beginner Aug 29 '18 at 14:31
0

I'd split the problem in two operations, as described in the comment to the question. This is possible since the first column is exactly the same for each file, and every file has the same number of lines.

Position yourself in the ballgown directory:

cd ballgown

As a first step, create an output file containing the first column:

cut -f6 TCGA-A7-A6VW-01A/t_data.ctab > out.tab

The bulk of the work is done by a combination of find and perl:

find -iname t_data.ctab -exec perl -i.bak -lane 'if($.==1){$ARGV=~/([-\w]+)\/.*$/;$f=$1} if(1..eof&&($n=$.)){$a[$.]=$F[11];$a[1]=$f;next}; print "$_\t$a[$.-$n]"' {} out.tab \;

Note: This is a destructive action; original files will be preserved with an added .bak extension.


Non-destructive version, making use of sponge (also, find was replaced by a for loop):

for F in */t_data.ctab; do perl -lane 'if(1..eof&&($n=$.)){$a[$.]=$F[11];$a[1]=$ARGV=~s/([-\w]+)\/.*$/$1/r;next} print "$_\t$a[$.-$n]"' $F out.tab | sponge out.tab; done;
0

Fully programmatic solution, in PHP.

<?php
$filenames = glob('*/t_data.ctab');
foreach($filenames as $k=>$filename) {
    $name = pathinfo($filename)['dirname'] . "\n";
    $file = file($filename);
    foreach ($file as $n => $line) {
        $line = explode("\t", $line);
        if ($n === 0) {
            $line[11] = $name;
        }
        if ($k === 0) {
            $out[$n] = $line[5] . "\t" . $line[11];
        } else {
            $out[$n] = trim($out[$n]) . "\t" . $line[11];
        }
    }
}
file_put_contents('out.tab', $out);

Usage:

  • Position yourself in the ballgown directory
  • Save the file with a name, let's say script.php
  • Run the script with php script.php
  • You'll find the output in the out.tab file

Note:

Let me know if you need further explanations regarding how to install and use PHP, what the script does and how to tweak it for a particular need.


Here's the same solution in Python, since the language was mentioned in the comments. This is the first time I write any Python, so please come forward with suggestions for improvements.

import os, glob
out = []
for k, filename in enumerate(glob.glob('*/t_data.ctab')):
    with open(filename, 'r') as f:
        file = f.readlines()
        for n, line in enumerate(file):
            line = line.split("\t")
            if n == 0:
                line[11] = os.path.dirname(filename) + "\n"
            if k == 0:
                out.append(line[5] + "\t" + line[11])
            else:
                out[n] = out[n].strip() + "\t" + line[11]
outfile = open('out.tab', 'w')
outfile.write("".join(out))

Same approach, written as a Perl one-liner:

perl -lane '$a[$n].=($a[$n]?"":$F[5])."\t".($n<1?$ARGV=~s#([-\w]+)\/.*$#$1#r:$F[11]); $n=eof?0:$n+1}{print "$_" for @a' */t_data.ctab
  • If i were to use scripting language, I would go with python. – zee Aug 29 '18 at 14:33
  • @zee That's a good idea, please provide a Python solution. – simlev Aug 29 '18 at 14:44

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.