Counting size of transcript from GTF file

Here I use grep to get records associated with query transcript and awk to perform the calculations:

cat test.gtf | grep "transcript_id \"ENSMUST00000105216\"" | egrep -v "CDS|start_codon|stop_codon" | awk '{l += $5 -$4 + 1}END{print l}'

Note that egrep is used to filter out all records containing CDS, start_codon or stop_codon. egrep is applied because it supports regular expression "logical or" statement.

Using regular expressions to analyze CIGAR

Suppose we want to know what is the value of a specific CIGAR operator in a SAM record. This problem can be easily solved by using regular expressions. I will use Python, but one can adapt code to any other programming language supporting RegExp.

For example, let's search for all possible alignments matches (M operator) in CIGAR:

In [55]: import re

In [56]: match = re.findall(r'(\d+)M', '40M25N5M')

In [57]: print match
-------> print(match)
['40', '5']


Expression \d+M represents all strings having pattern "nM", where n is a number that can consist from one or multiple digits. Round brackets create a group from the number, so it can be accessed later.

Similarly one can iterate over a CIGAR string:


In [60]: match = re.findall(r'(\d+)(\w)', '40M25N5M')

In [61]: match
Out[61]: [('40', 'M'), ('25', 'N'), ('5', 'M')]


Here we use \w meta-symbol to represent any letter and round brackets for grouping.

Have fun!