1

I have a large GTF file, like below:

 # ./stringtie -p 4 -G /home/humangenome_hg19/homo_gtf_file.gtf -o strAD1_as/transcripts.gtf -l strAD1 /home/software/star-2.5.2b/bin/Linux_x86_64/mapA1Aligned.sortedByCoord.out.bam                               
# StringTie version 1.3.2d                              
1   StringTie   transcript  30267   31109   1000    +   .   gene_id "strAD1.1"; transcript_id "strAD1.1.1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.028725"; FPKM "0.053510"; TPM "0.109957";
1   StringTie   exon    30267   30667   1000    +   .   gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.014218";
1   StringTie   exon    30976   31109   1000    +   .   gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "2"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.072139";

I want to have the 9th column with just gene_id, transcript_id, reference_id and ref_gene_id. They are in the 9th column and separated by space (the columns themselves are TAB-separated). Could you please help me out how I can such a column with a simple command in Linux? I don't want to use Excel for it.

2
  • The rest of the column are delimited by tabs?
    – Kusalananda
    May 13, 2017 at 7:40
  • Yes, they are delimited by tabs
    – Mary
    May 13, 2017 at 7:46

2 Answers 2

3

Ideally, since the data is in GTF format, one should use a GTF parser to parse it. I currently have no such parser or parsing library installed so my solution is based solely on the data that you have provided in the question.

To extract the 9th column:

$ cut -f 9 data.gtf
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.028725"; FPKM "0.053510"; TPM "0.109957";
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "1"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.014218";
gene_id "strAD1.1"; transcript_id "strAD1.1.1"; exon_number "2"; reference_id "ENST00000469289"; ref_gene_id "ENSG00000243485"; ref_gene_name "MIR1302-10"; cov "0.072139";

To get the data that we want from this, we need to treat transcripts and exons separately as their attributes have different order in the data. We do this with awk and output different fields in the input data depending on whether the current line contains the string exon_number or not:

$ cut -f 9 data.gtf | awk '/exon_number/ { print $2, $4, $8, $10; next } { print $2, $4, $6, $8 }'
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";
"strAD1.1"; "strAD1.1.1"; "ENST00000469289"; "ENSG00000243485";

Then we remove the double quotes and semicolons from this:

$ cut -f 9 data.gtf | awk '/exon_number/ { print $2, $4, $8, $10; next } { print $2, $4, $6, $8 }' | tr -d '";'
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
strAD1.1 strAD1.1.1 ENST00000469289 ENSG00000243485
3

Maybe just:

< file cut -sd '"' -f2,4,8,10 | tr '"' ' '

That is consider the input as a list of "-separated columns and extracting the 2nd, 4th, 8th and 10th columns.

With GNU cut, you can replace the | tr '"' ' ' with --output-delimiter=' '.

That makes the assumption that " characters don't appear elsewhere in the lines, that those gene_id, transcript_id... attributes always appear and always in that order.

As noted by Kusalananda, that's not the case in your sample where it should be 2,4,6,8 for the first line and 2,4,8,10 for the others.

To do a more expressive matching: that only the 9th tab-delimited column should be considered and that the correct attribute names be found, you could resort to regular expressions like with:

< file pcregrep -o1 -o2 -o3 -o4 --om-separator=' ' '(?x)
  ^(?:[^\t]*+\t){8}(?=[^\t]*? \b gene_id       \ +"([^"\t]*)")
                   (?=[^\t]*? \b transcript_id \ +"([^"\t]*)")
                   (?=[^\t]*? \b reference_id  \ +"([^"\t]*)")
                   (?=[^\t]*? \b ref_gene_id   \ +"([^"\t]*)")'

If you don't have pcregrep or too old a version to support -o1..., you can use perl instead:

< file perl -lne 'print "$1 $2 $3 $4" if m{
  ^(?:[^\t]*+\t){8}(?=[^\t]*? \b gene_id       \ +"([^"\t]*)")
                   (?=[^\t]*? \b transcript_id \ +"([^"\t]*)")
                   (?=[^\t]*? \b reference_id  \ +"([^"\t]*)")
                   (?=[^\t]*? \b ref_gene_id   \ +"([^"\t]*)")}x'

That regexp first matches the first 8 fields ((?:[^\t]*+\t){8}) and following that, we have 4 look-ahead expressions ((?=...)), so we're matching for those 8 fields provided what follows matches all 4 of the look-ahead expressions. Each look-ahead expression looks for one of the attributes and captures the value (in the (...) part). Those captured values are then available in $1, $2, $3, $4.

So that allows attributes in any order.

Note that it could be fooled by things like:

1 2 3 4 5 6 7 8 gene_id "transcript_id " ...

While it would be possible to address it, it's probably not worth the effort as I don't expect it would be occurring in the input.

While you're at using perl, you could also do a more formal parsing of that 9th field. Something like:

< file perl -F'\t' -lane '
  my %field;
  while ($F[8] =~ /(\w+) +"(.*?)"/g) {$field{$1}=$2}
  if (%field) {
    print join " ", @field{
      qw(gene_id transcript_id reference_id ref_gene_id
    )}
  }'

(here, printing a line as long as at least one attribute is found (as opposed to all the requested attributes in the other approaches)).

1
  • Thanks for all responses. Both commands worked fine. but I prefer the first one.
    – Mary
    May 13, 2017 at 10:20

You must log in to answer this question.

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