-1

I have 6 text files (each corresponds to a specific sample) and each file looks like this:

Gene_ID Gene_Name   Strand  Start   End Length  Coverage    FPKM    TPM
ENSMUSG00000102735  Gm7369  +   4610471 4611406 936 0   0   0
ENSMUSG00000025900  Rp1 -   4290846 4409241 10926   0   0   0
ENSMUSG00000104123  Gm37483 -   4363346 4364829 1484    0   0   0
ENSMUSG00000102175  Gm6119  -   4692219 4693424 1206    0.328358    0.015815    0.008621

I want to collect all the elements from 1 & 2 column in one file and corresponding tpm values(9th column) for each sample in a new file, so wherever there is no tpm value enter 0.

My output file should look like this:

gene_id gene_name sample1_tpm sample2_tpm sample3_tpm ......sample6_tpm
0

This solution uses awk to gather TPM values from all 6 input files file1file6 in an array indexed by gene_id and gene_name, and finally prints the consolidated TPM values for each array index. Missing TPM values are replaced by 0.

awk ' BEGIN { f=0 } { if ($1 == "Gene_ID") { f++; getline; } i = $1 " " $2 row[i] = row[i] " " (($9 == "") ? "0" : $9) } END { for (i in row) print i,row[i] } ' file1 file2 file3 file4 file5 file6

0

Here's a perl script that does what you want. It builds up a hashed array (indexed by a string made up of "gene_id gene_name") containing an array of tpm samples.

Once it's finished reading all the files, it prints out each element of the hashed array. Note the script does not require exactly 6 elements, nor is it limited to <= 6 - it will work for any number of tpm elements.

If some gene_id's have 6 tpm array elements while others have 1 or 5 or 1000 elements, each gene_id will be printed with the exact number of tmp elements detected for that gene_id/gene_name. See below for extra code you can add to treat this as an error rather than a feature.

Input can be provided on stdin and/or by giving it filenames to process on the command line.

#! /usr/bin/perl

my %samples = () ;

while(<>) {
    chomp;
    s/^\s*|\s*$//g;
    next if (m/^Gene_ID|^$/);

    my ($gene_id, $gene_name, undef,undef,undef,undef,undef,undef, $tpm) = split;
    $tpm = 0 unless (defined($tpm));

    push @{ $samples{"$gene_id $gene_name"} }, $tpm;
}

foreach my $key ( keys %samples ) {
    my ($gene_id, $gene_name) = split(/\s+/,$key);

    print "$gene_id\t$gene_name\t", join("\t", @{ $samples{"$gene_id $gene_name"} }), "\n";
}

The script's output is tab-delimited. If you prefer something else, edit the print line in the script and replace ALL the \t's with your preferred delimiter.

If you want to produce an error message and exit if there are not exactly 6 tpm elements printed for each gene_id, add the following code immediately after the my ($gene_id, $gene_name) = split(/\s+/,$key) line:

    my $count = scalar @{ $samples{"$gene_id $gene_name"} };
    if ($count != 6) { die "Warning: $gene_id $gene_name has $count element(s)\n"};

Sample Output:

I made 6 copies of the sample data you provided, called sample.txt to sample6.txt. Saved the perl script above as sample.pl, made it executable with chmod +x sample.pl, and ran it:

$ ./sample.pl sample*.txt
ENSMUSG00000104123  Gm37483 0   0   0   0   0   0
ENSMUSG00000102175  Gm6119  0.008621    0.008621    0.008621    0.008621    0.008621    0.008621
ENSMUSG00000102735  Gm7369  0   0   0   0   0   0
ENSMUSG00000025900  Rp1 0   0   0   0   0   0

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.