Spring Into Bioinformatics

Gene Annotation

In the last session we looked at our aligned sample RNA (cDNA) fragments to the Mouse GRCm38/mm10 reference genome which contains the location of each fragment and the quality of each alignment. However the alignment process will only tell you where a RNA fragment has mapped, but not what it has mapped too. The primary of an RNA-seq experiment is to quantify the level of expression of each gene, so in this session we are going to use the reference gene annotation to add context to our previous alignments.

If the genome is a road in which you drive on, then the genome annotation would be the map that identifies where everything is located. For example, the Foxp3 (forkhead box P3) gene, an important transcription factor in the regulation of T cells in mice, is located on the X chromosome between 7,579,676bp and 7,595,245bp. We can use this information to count all the reads that align to that region to establish the quantity of RNA we have in our sample.

Genome annotations not only contain information about genes, but all other genomic features across the genome. Locations of things like transposable elements (TEs) and repeats, centromere and telomere locations and non-coding RNAs can all be contained within these files

GFF/GTF Files

There are many files that can contain gene annotation information however today we will be using General Feature Format (GFF) or General Transfer Format (GTF) files. GFF files have version 2 and version 3 formats, which are slightly different and today, we’ll just look at GTF files, which are best considered as GFF2.2. The differ only by the restrictions that are placed on the type of entries that can be placed in some columns.

GTF files are basically text files with columns that contain standard information. Columns are tab-separated with no line provided that gives the column names. These are fixed by design, and as such, explicit column names are not required.

1. seqname - name of the chromosome or scaffold; chromosome names can be given with or without the ‘chr’ prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
2. source - name of the program that generated this feature, or the data source (database or project name)
3. feature - feature type name, can only take the values “CDS”, “start_codon”, “stop_codon”, “5UTR”, “3UTR”, “inter”, “inter_CNS”, “intron_CNS” and “exon” (CNS stands for Conserved Noncoding Sequence)
4. start - Start position of the feature, with sequence numbering starting at 1.
5. end - End position of the feature, with sequence numbering starting at 1.
6. score - A floating point value (i.e. decimal points are allowed)
7. strand - defined as + (forward) or - (reverse).
8. frame - One of ‘0’, ‘1’ or ‘2’. ‘0’ indicates that the first base of the feature is the first base of a codon, ‘1’ that the second base is the first base of a codon, and so on…
9. attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature. In the GTF format two mandatory features are required here, although they can be an empty string, i.e. "":
   • gene_id value
   • transcript_id value

An example is given below. Note that header rows are not controlled, but must start with the comment character #

# Data taken from http://mblab.wustl.edu/GTF22.html
381 Twinscan  CDS          380   401   .   +   0  gene_id "001"; transcript_id "001.1";
381 Twinscan  CDS          501   650   .   +   2  gene_id "001"; transcript_id "001.1";
381 Twinscan  CDS          700   707   .   +   2  gene_id "001"; transcript_id "001.1";
381 Twinscan  start_codon  380   382   .   +   0  gene_id "001"; transcript_id "001.1";
381 Twinscan  stop_codon   708   710   .   +   0  gene_id "001"; transcript_id "001.1";

Note: People variously use GFF and GTF to talk about GFF version 2, and GFF to talk about GFF version 3. GFF2 is not compatible with GFF3, so make sure you have the correct file format if you are given a GFF file. There are conversion tools available to inter-convert them, they are rarely reliable.

Mouse Chromosome 1

For today, we’ll need the GTF from MouseChromosome 1. Tomorrow we will also need the Mouse Chromosome 1 sequence, so lets download them all and put them in a location on our VM. We’ll get these from Ensembl, so cut and paste the following code directly into your VM.

Note: For those of you that are unfamilar with coding, in bash we can add comments to code to explain what we are doing. This is good practise for anyone working on the command-line, as its open easy to forget why you are running something. Any line starting with a # will be ignored by the program.

cd

# Create a directory for the mouse genome and change into it
mkdir -p genomes/Mmusculus
cd genomes/Mmusculus

# Dowload the Chromosome 1 fasta file from a weblink
wget ftp://ftp.ensembl.org/pub/release-93/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.1.fa.gz

# Uncompress the file (its gzipped at the moment) so we can use it later
gunzip Mus_musculus.GRCm38.dna.chromosome.1.fa.gz

# Download the genome annotation (but keep this gzipped up)
wget ftp://ftp.ensembl.org/pub/release-93/gtf/mus_musculus/Mus_musculus.GRCm38.93.chr.gtf.gz

You may have already realised that this gtf contains all the information for the entire genome, rather than just chromosome 1, so we’ll need to create one containing just the gene information for that sequence.

First we’ll copy the file header.

zcat Mus_musculus.GRCm38.93.chr.gtf.gz | egrep "^#!" > Mus_musculus.GRCm38.93.chr1.gtf

Then we’ll just get the gene information for chr1. The first field in a gtf file is the chromosome, we this will be easy to extract using egrep

zcat Mus_musculus.GRCm38.93.chr.gtf.gz | egrep "^1\s" >> Mus_musculus.GRCm38.93.chr1.gtf

To avoid confusing ourselves, let’s delete the full gtf so we only have the one.

rm Mus_musculus.GRCm38.93.chr.gtf.gz

Questions

1. How many genes and exons are found on chromosome 1?

1. Declare variables

Here we need to declare variables that will be used in the command, such as the number of cores that the program will run parallel, the file path of the GTF file that we downloaded above and the directory that contains our alignment data. Generally we run this as a script on the command-line, however today we have just run it on the terminal interactively.

2. Find our samples

Note that we’ve used a different strategy for finding our files this time. Instead of using $(ls ...), we’ve used find. This is probably a superior approach and can withstand difficult file paths more easily. Output from find will give you each result on a new line, so after this we’ve used tr to translate line breaks \n to spaces.

3. Count our alignments

Here we’ve called featureCounts and passed several parameters to it. The help page is quite extensive for this tool, so feel free to browse it in your terminal or on pages 37-42 of the manual pdf.

4. Tidy up our output

The output from featureCounts contains a few columns which we ignore, even though others may find them useful (exon coordinates, gene lengths etc) In the final line, we’re using cut to just return columns 1, then everything beyond 7. sed is then used to delete the first line. This will give us data in a nice easy format for importing into R.

Tasks

• Run featureCounts with the an unstranded library parameter (-s 1). How do they differ?