Introduction to Unix for Bioinformatics

Handling genomic coordinates

Contents

• Opening a session
• Downloading RefSeq transcript coordinates
• Some descriptive statistics on refSeq transcript using basic unix command
• The Bedtools suite
• Getting promoter regions of each transcript
• Intersecting promoter regions with predicted transcription factor binding sites
• Functional annotation using Database for Annotation, Visualization and Integrated Discovery

Opening a session.

Starting a terminal

Using this terminal you should now be able to communicate with the system through a set of shell commands (your default shell command interpreter should be BASH).

Creating a working directory

	cd ~             # or cd
	pwd              # check the current working directory
	mkdir td1_genomic_coordinates
	cd td1_genomic_coordinates # or cd ~/td1_genomic_coordinates
      

Downloading RefSeq transcript coordinates

RefSeq transcript

The Reference Sequence (RefSeq) collection aims to provide a comprehensive, integrated, non-redundant, well-annotated set of sequences. Here we will download informations about RefSeq transcripts (that is known full-length transcript) related to Mus musculus genome (mm9 version). Refseq transcript can be obtained from UCSC genome browser.

	
	wget http://hgdownload.cse.ucsc.edu/goldenPath/mm9/database/refGene.txt.gz
	gunzip refGene.txt.gz
      

The UCSC Web page describes the column content of the refGene table that we just downloaded. This table contains one line per "object" (gene), and one column per attribute. This description table is of course crucial to understand the content of the refGene.txt file that we downloaded in the previous step.

Let's now have a look at the refGene file.

		   
	wc -l refGene.txt
	head -n 5 refGene.txt
	tail -n 5 refGene.txt
	less refGene.txt
	head -n 1 refGene.txt  | od -c 					
      

Some descriptive statistics on refSeq transcript using basic unix command

	cut -f13 refGene.txt
	cut -f13 refGene.txt | grep "^Bcl"
	cut -f13 refGene.txt | grep "^Bcl" | sort | uniq | wc -l		
	cut -f13 refGene.txt | sort | uniq -c | sort -n | tail -1
	cut -f13 refGene.txt | sort | uniq -c | sort -n | tail -1
	cut -f9,13 refGene.txt | sort -n | tail -1
	cut -f2 refGene.txt  | sed 's/_.*//' | sort | uniq -c
	cut -f3 refGene.txt  | sort | uniq -c | sort -n
      

The Bedtools suite

What is bedtools

The bedtools manual says:

The BEDTools utilities allow one to address common genomics tasks such as finding feature overlaps and computing coverage. The utilities are largely based on four widely-used file formats: BED, GFF/GTF, VCF, and SAM/BAM. Using BEDTools, one can develop sophisticated pipelines that answer complicated research questions by "streaming" several BEDTools together. The following are examples of common questions that one can address with BEDTools.

1. Intersecting two BED files in search of overlapping features.
2. Computing coverage for BAM alignments based on genome features.
3. Merging overlapping features.
4. Screening for paired-end (PE) overlaps between PE sequences and existing genomic features.
5. Calculating the depth and breadth of sequence coverage across defined "windows" in a genome.
6. Screening for overlaps between "split" alignments and genomic features.

Installing the bedtools suite

Check if the bedtools suite is already installed on your system (type intersectBed in a terminal). If bedtools is not available, follow the procedure below to perform installation.

      wget https://bedtools.googlecode.com/files/BEDTools.v2.16.2.tar.gz
      mkdir -p  ~/bin
      mv BEDTools.v2.16.2.tar.gz ~/bin
      cd ~/bin
      tar xvfz BEDTools.v2.16.2.tar.gz
      cd BEDTools-Version-2.16.2
      make
      echo PATH=$PATH:~/bin/BEDTools-Version-2.16.2/bin >> ~/.bashrc
      source ~/.bashrc
    

Transforming refGene file into bed format

First we need to transform the refGene.txt file into a bed to use it as input for the bedtools suite. We will do it with awk.

      awk 'BEGIN{FS=OFS="\t"}{print $3,$5,$6,$2"|"$13,0,$4}' refGene.txt> refGene.bed
    

Upload the resulting bed file into UCSC (top menu, genome > add custom track > choose file > submit). Go to the genome browser to visualize the new track.

Getting promoter regions of each transcript

To obtain promoter regions of each transcript we will first extract the TSS coordinates of each transcript. The TSS will correspond to the txStart or txEnd coordinates depending on wether genes is located on the positive or negative strand respectively ( see refSeq table definition).

      awk 'BEGIN{FS=OFS="\t"}($6=="+"){print $1,$2,$2,$4,$5,$6}' refGene.bed > refGene.TSS.bed
      awk 'BEGIN{FS=OFS="\t"}($6=="-"){print $1,$3,$3,$4,$5,$6}' refGene.bed >> refGene.TSS.bed
    

Now we will use slopBed to increase the size of each feature (TSS) to define a promoter region around the TSS ([-2000,+500]). Also this can be done with a simple awk command, sploBed will restrict the resizing to the size of the chromosome (i.e. no start < 0 and no end > chromosome size). We thus need to download the length of each chromosome (here from UCSC ftp site).

      wget http://hgdownload.cse.ucsc.edu/goldenPath/mm9/database/chromInfo.txt.gz
      gunzip chromInfo.txt.gz
      cat chromInfo.txt
    

Now we can run slopBed. We will also reorder features (TSS) according to chromosomes coordinates

      slopBed -i refGene.TSS.bed -l 2000 -r 500 -s  -g chromInfo.txt  | sortBed > refGene_Prom_2000_500.bed
    

Check the results by loading the refGene_Prom_2000_500.bed as a new user track in the UCSC genome browser.

NB: The RSAT web server also proposes advance tools for retrieving promoter regions (excluding coding regions, repeat elements,...). This tools can also be run using UNIX command lines.

Intersecting promoter regions with predicted transcription factor binding sites

The dataset that will be used thereafter was obtained from the following publication:

• Lepoivre C, Bergon A, Lopez F, Perumal NB, Nguyen C, Imbert J, Puthier D. TranscriptomeBrowser 3.0: introducing a new compendium of molecular interactions and a new visualization tool for the study of gene regulatory networks. BMC Bioinformatics. 2012 Jan 31;13:19.

In this article, the authors scanned vertevrates alignments in search for highly conserved putative transcription factor binding sites (TFBS). The resulting file contains the coordinates of the predicted TFBS (mm9 genome version)

	wget http://www.bigre.ulb.ac.be/courses/statistics_bioinformatics/data/cis-reg/conserved_predicted_TFBS_mm9.bed
	head conserved_predicted_TFBS_mm9.bed
	intersectBed -wa -wb -a refGene_Prom_2000_500.bed -b conserved_predicted_TFBS_mm9.bed > conservedsites_prom_intersect.bed
	grep  -w M00920 conservedsites_prom_intersect.bed   | sed 's/.*|//' | cut -f1 | sort | uniq > M00920_targets.txt
	grep  -w M00771  conservedsites_prom_intersect.bed  | sed 's/.*|//' | cut -f1 | sort | uniq > M00771_targets.txt
	grep  -w M00941  conservedsites_prom_intersect.bed  | sed 's/.*|//' | cut -f1 | sort | uniq >  M00941_targets.txt
	grep  -w M00054  conservedsites_prom_intersect.bed  | sed 's/.*|//' | cut -f1 | sort | uniq >  M00054_targets.txt
	grep  -w M01028  conservedsites_prom_intersect.bed  | sed 's/.*|//' | cut -f1 | sort | uniq >  M01028_targets.txt
      

Functional annotation using Database for Annotation, Visualization and Integrated Discovery

We have seen that each transcription factor is associated with a set of putative targets. What can we learn from these targets? Is the transcription factor associated to genes involved in some particular biological functions?

To answer this question we will use the Database for Annotation, Visualization and Integrated Discovery (DAVID). DAVID will try to find whether the selected genes share some particular features, by analyzing their annotations from several sources including Gene Ontology, KEGG pathways, Reactome Pathways, genomic locations, ...

Note: the Fisher test will be explained in the next practical "GO statistics".