Skip to content

Instantly share code, notes, and snippets.

@darencard
Last active March 28, 2024 13:30
Show Gist options
  • Save darencard/9497e151882c3ff366335040e20b6714 to your computer and use it in GitHub Desktop.
Save darencard/9497e151882c3ff366335040e20b6714 to your computer and use it in GitHub Desktop.
Extracting spliced sequences (e.g., CDS) from GFF files

Extracting spliced sequences (e.g., CDS) from GFF files

GFF is a common format for storing genetic feature annotations. In the case of gene annotations, subsets of elements are split over multiple lines, as things like exons and CDS features will have gaps based on the full genome sequence. Therefore, while it is easy to extract exon and CDS lines, it can be difficult to associate them together based on a parent (e.g., transcript) ID and perform downstream operations. Even extracting the full CDS sequence using a GFF file can be tricky for this reason, even though it seems trivial.

Here we'll overcome this difficulty using the gffread tool. Installation is pretty easy and is documented in the GitHub README. gffread has a lot of options, but here we'll just document one that extracts the spliced CDS for each GFF transcript (-x option). Note that you can do the same thing for exons (-w option) and can also produce the protein sequence (-y option).

Let's extract the CDS sequences for each transcript using a genome sequence and a GFF annotation file.

gffread -x <out.fasta> -g <genome.fasta> <annotation.gff>

Pretty straight-forward. We can also set the command to output to STDOUT instead of a file by seeing the option -x -. Let's demo this by going a step further with the output and extracting every 3rd codon position (relies on bioawk). Key is an awk command documented here.

gffread -x - -g <genome.fasta> <annotation.gff> | \
bioawk -c fastx '{ print $name, $seq }' | \
while read line; \
do \
name=$(echo $line | cut -f 1); \
echo $line | cut -f 2 | \
awk -F "" '{ for (i = 3; i <= NF; i += 3) \
printf "%s%s", $i, (i+3>NF?"\n":FS) }' | \
awk -v name="$name" '{ print ">"name; print $1 }'; \
done \
> <out.fasta>

Now it is easy to calculate something like GC3, the GC content of 3rd codon positions, using a tool like seqtk. We'll make the output in BED format.

gffread -x - -g <genome.fasta> <annotation.gff> | \
bioawk -c fastx '{ print $name, $seq }' | \
while read line; \
do \
name=$(echo $line | cut -f 1); \
echo $line | cut -f 2 | \
awk -F "" '{ for (i = 3; i <= NF; i += 3) \
printf "%s%s", $i, (i+3>NF?"\n":FS) }' | \
awk -v name="$name" '{ print ">"name; print $1 }'; \
done | \
seqtk comp - | \
awk -v OFS="\t" '{ print $1, "0", $2, ($4 + $5) / $2 }'

It is even easier if you want to look at GC across the entire CDS.

gffread -x - -g <genome.fasta> <annotation.gff> | \
seqtk comp - | \
awk -v OFS="\t" '{ print $1, "0", $2, ($4 + $5) / $2 }'

Analytics

@thibui1996
Copy link

Hello,

In the section where you said you would make the output in BED format, where in the code did you indicate the output to be bed? Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment