Exercise 2. Function enrichment analysis & gene function annotation

After RNA-seq data analysis, you will get a list of differentially expressed (DE) genes. In this exercise, we will use two different ways to test which function categories are enriched in your DE gene list.

If you work on a non-model organism, chances are that you do not have Gene Ontology (GO) annotation of your reference genome. GO annotation is required for function enrichment analysis. We will use BLAST2GO software to generate GO annotation.

There are three GO domains: cellular component (CC), biological process (BP), and molecular function (MF). In this exercise, we will only do function enrichment test for the BP domain. When you do your real project, you will need to test all three.

Part 1. Over representation analysis (ORA)

The topGO package of Bioconductor will be used to do ORA. We will perform ORA on RNA-seq results of a yeast experiment. The reference is https://www.ncbi.nlm.nih.gov/pubmed/16606683 .

  1. Download GO annotation from Ensembl BioMart.

    Ensembl (https://ensembl.org/) is a good resource to retrieve existing GO annotation (For plant genomes, go to https://plants.ensembl.org/ )

  2. Convert the file you downloaded into a format compatible with topGO. Use FileZilla to upload the "mart_export.txt" file to your assigned Linux server, keep it in the /workdir/$USER.

We provide you with a script that converts the BioMart file into topGO format. Run the commands:

cd /workdir/$USER
/shared_data/annotation2019_2/toTopGO.pl mart_export.txt

After this, you would see a new file created: topgoAnnot.

  1. You are provided with two gene lists: refset and DElist_pos. Run the R script topGO.r:
   cp /shared_data/annotation2019_2/* ./
   Rscript /shared_data/annotation2019_2/topGO.r topgoAnnot refset DElist_pos 0.1 BP myBP

· topgoAnnot: the GO annotation file in the topGO required format;

· refset: a text file with list of reference set of genes with one gene per line (normally all genes that have none-zero expression in your experiments)

· DElist_pos: a text file with the list of up-regulated DE genes.

· 0.1: cutoff p-value for enriched categories.

· BP: test for biological processing GO. You can also test for MF (molecular function) and CC (cellular component).

· myBP: output file name prefix.

After this, you should find two new files "myBP".

You can use FileZilla to download both files to your laptop to examine the contents.

In the myBP file, the enriched GO ids are listed sorted by p-values in the column "topgoFisher".

In your publication, you can say that "gene set enrichment analysis was performed with BioConductor topGO package, using its default method weight01, and p-value threshold is set at ...".


Part 2. Gene Set Enrichment Analysis (GSEA)

We provide you with two data files (bp.gmt & genes.rnk) to run GSEA. When you work on your own project, here is how to prepare the two files.

You can either GSEA on your laptop, or on the LINUX server.

  1. Install and run GSEA on your laptop.
  1. Examine results Click "Success". If it does not work, click "Show results folder", which should open the directory. Open the directory of your analysis, double click the file "index.html".

On the page, you should see two block of "Enrichment in phenotype". One for gene sets enriched in up-regulated genes (na_pos) and one for gene sets enriched in down-regulated genes (na_neg). Click "Detailed enrichment results in html ", you would see enriched gene sets.

  1. (Optional) Enrichment map visualization. Interpretation of the map can be found at this page: https://enrichmentmap.readthedocs.io/en/docs-2.2/Tutorial_GSEAInterface.html

  2. Run GSEA command on LINUX (optional).

With GSEA open on your laptop, click "command" at the bottom of the GSEA window. Copy-paste the command to a text editor.

You need to change a few things for this command to run on BioHPC computer:

Now run the command on Linux server. After done, you can examine the results in the "out" directory.

Part 3. Run BLAST2GO to create GO annotation. (Do it at home)

In part 1 of this exercise, you have copied all files in /shared_data/annotation2019_2/ to /workdir/$USER. You will run BLAST2GO on the fasta file: annot_exercise_aa.fasta

  1. Run BLASTX against SWISSPROT database. BLASTX is a command line software.
  cd /workdir/$USER
  cp /shared_data/genome_db/BLAST_NCBI/swissprot* ./
  blastp -num_threads 4 -query annot_exercise_aa.fasta  -db swissprot \
    -out blastresults.xml  -max_target_seqs 20 -evalue 1e-5 -outfmt 5 \
    -culling_limit 10 >& logfile &
  1. After it is done, copy the result file into your home directory.
cp blastresults.xml   ~/
  1. Run BLAST2GO on cbsumm10.

    As we only have one license of BLAST2GO. It is installed on the computer cbsumm10.biohpc.cornell.edu.

    Login to the computer cbsumm10.

mkdir /workdir/$USER
cd /workdir/$USER
cp /shared_data/blast2go/* ./
cp ~/blastresults.xml  ./
/usr/local/blast2go/blast2go_cli.run -properties annotation.prop -useobo go.obo \
     -loadblast blastresults.xml -mapping -annotation -annex -statistics all \
     -saveb2g myresult -saveannot myresult -savereport myresult -tempfolder ./ \ 
     >& annotatelogfile &

After this step, you will get three result files: