Jetstream Workflow That Creates Gene Expression Matrices (GEMs) from SRA/FASTQ NGS Files

Rationale and background

Pegasus-GEM is Pegasus Jetstream workflow that utilizes Jetstream resources to produce a Gene Expression Matrix (GEM) from DNA sequence files in FASTQ format. This is adapted from OSG-GEM workflow that currently runs on Open Science Grid (OSG). 

William L. Poehlman, Mats Rynge, Chris Branton, D. Balamurugan and Frank A. Feltus. OSG-GEM: Gene Expression Matrix Construction Using the Open Science Grid. Bioinformatics and Biology Insights 2016:10 133–141 doi: 10.4137/BBI.S38193.

Introduction

This workflow processes both paired or single end FASTQ files to produce a matrix of normalized RNA molecule counts (FPKM). Pegasus-GEM also supports direct input downloads from NCBI SRA for processing. An indexed reference genome along with gene model annotation files must be obtained prior to configuring and running the workflow. The following tasks are directed by the Pegasus workflow manager:

• Splitting input FASTQ files into files containing 20,000 sequences each.
• Trimming raw sequences with Trimmomatic
• Aligning reads to the reference genome using Hisat2 or Tophat2 or STAR
• Merging BAM alignment files into a single sorted BAM file using Samtools merge
• Quantifying RNA transcript levels using StringTie or Cufflinks

It is suggested that the user become familiar with the documentation associated with the following software packages:

• Trimmomatic
• Hisat2
• Tophat2
• STAR
• Bowtie2
• Samtools
• StringTie
• Cufflinks

This tutorial will take users through steps of:

1. Running Pegasus-GEM on Jetstream cloud
2. Running Pegasus-GEM on an example data

Considerations

Sounds great, what do I need to get started?

1. XSEDE account
2. XSEDE allocation. Users can request a startup XSEDE allocation. The startup XSEDE allocations are relatively easy to apply and get allotted. Please contact help@xsede.org for further information.
3. Getting started with Jetstream
4. Your data (or you can run example data)

What kind of data do I need?

1. Mandatory requirements
   1. FASTQ files. Either single end or paired. You can also use SRR id's
   2. Select the type of reads. Single or Paired.
   3. Select the type of mapper (Tophat2 or Hisat2 or STAR)
   4. Select the type of assembler (Cufflinks or Stringtie)

What kind of resources will I need for my project?

1. Enough SUs (Service Units) to run your computation
2. One MASTER and several workers needed for running your computation
3. Enough storage space on the WQ-MAKER Jetstream instance for both input and output files. Since most of the images have limited disk space, it is recommended to mount an external volume to the running Pegasus-GEM MASTER instance. 

Part 1: Connect to an instance of an CentOs Jetstream Image (virtual machine)

Step 1. Go to https://use.jetstream-cloud.org/application and log in with your XSEDE credentials.

Step 2. Click on the "Create New Project"  in the Project tab on the top and enter the name of the project and a brief description 

Step 3. Launch an instance from the selected image and name it as MASTER 

After the project has been created and entered inside it, click the "New" button, select "Centos 7 (7.3) Development GUI" image and then click Launch instance. In the next window (Basic Info),

• name the instance as "MASTER" (don't worry if you forgot to name the instance at that point, as you can always modify the name of the instance later)
• leave base image version as it is
• leave the project name as it is or change to a different project if needed
• select "Jetstream - Indiana University or Jetstream - TACC" as Provider (for this tutorial we will chose Jetstream - TACC) and click 'Continue'. Your choice of provider will depend on the resources you have available (SUs) and the needs of your instance
• select "m1.medium" as Instance size and click "Continue". 
   

Launch one instance from the  "Centos 7 (7.3) Development GUI" image and name them as W-1. Use the same configuration as the MASTER VM 
 

$ ssh -A root@<ip.address>
# For example
$ ssh -A root@129.114.17.50 #Do not use this ip address with your username as it will not work

[js-17-50] root ~-->
cd /srv/
git clone https://github.com/pegasus-isi/jetstream-pegasus.git
cd jetstream-pegasus/
./bootstrap-master.sh # MASTER     
./bootstrap-worker.sh 129.114.17.34 # W-1
exit

Since the m1 medium instance size (60GB disk space) selected for running MASTER instance of Pegasus may not be sufficient for most of the GEM runs, it is recommended to run it on volumes. Following are the three steps for mounting the external volume.

6.1: Create a volume

Click the "New" button in the project and select "Create Volume". Enter the name of the volume, volume size (GB) needed and the provider (TACC or Indiana) and finally click "Create Volume" 

6.2: Attach the created volume to the MASTER instance 

Click the checkbox by the side of Volume (1), then click the attach button (2), next chose the MASTER instance (3) and finally click the ATTACH VOLUME TO INSTANCE button (4)

# Create a directory .ssh if you don't have that already
$ mkdir -p ~/.ssh 
$ ssh-keygen -t rsa -b 2048 -f ~/.ssh/workflow 
(Just hit enter when asked for passphrase)

Generating public/private rsa key pair.
Enter passphrase (empty for no passphrase): 
Enter same passphrase again: 
Your identification has been saved in /home/upendra/.ssh/workflow.
Your public key has been saved in /home/upendra/.ssh/workflow.pub.
The key fingerprint is:
e4:63:31:61:46:ef:52:0e:dd:04:70:bc:a8:46:0a:51 upendra@js-17-50.jetstream-cloud.org
The key's randomart image is:
+--[ RSA 2048]----+
|   .E  .*oo..    |
|  .    o =.o     |
|   .    =.+..    |
|  .   .o.B.      |
|   . o .S o      |
|    . o. o       |
|     .           |
|                 |
|                 |
+-----------------+
$ cat ~/.ssh/workflow.pub >>~/.ssh/authorized_keys

Step 2: Example Workflow Setup. The workflow cloned from github contains an example config file as well as example input files from the 21st chromosome of Gencode Release 24 of the GRCh38 build of the human reference genome. Two small FASTQ files containing 200,000 sequences from NCBI dataset SRR1825962 lie within the Test_data directory of the workflow

$ git clone https://github.com/feltus/OSG-GEM.git

$ cd OSG-GEM

# To run the test workflow, the user must copy the osg-gem.confi.template file like this
$ cp osg-gem.conf.example osg-gem.conf

# Take a look at the config file first
$ cat osg-gem.conf
[reference]
reference_prefix = chr21-GRCh38

[inputs]
#
# List the inputs to process. Each line can either be a pair
# of forward and reverse files, separated by a space:
#
#    input1 = forward.fastq.gz reverse.fastq.gz
#
# or a single SRR number. Example:
#
#    input2 = DRR046893
#
# or a single fastq file for single end reads.  Example:
#    input3 = SRR4343300.fastq.gz

input1 = ./Test_data/TEST_1.fastq.gz ./Test_data/TEST_2.fastq.gz
#input1 = ./Test_data/SRR4343300.fastq.gz
#input2 = SRR4343300
#input2 = DRR046893

[config]

# Memory available to the jobs. This should be roughly 2X the
# size of the reference genome, rounded up whole GB
memory = 4 GB

# Reads are single end
single = False

# Reads are paired end
paired = True

# process using TopHat2
tophat2 = False

# process using Hisat2
hisat2 = True

# process using STAR
star = False

# process using Cufflinks
cufflinks = False

# process using StringTie
stringtie = True

Run default test data . The workflow, configured to run Hisat2 and StringTie, can then be launched by running:

$ ./submit

# To check the status of your run. For example, run
$ pegasus-status -l /mydata/OSG-GEM/runs/osg-gem-1498085153/workflow/osg-gem-1498085153 
(no matching jobs found in Condor Q)
UNRDY READY   PRE  IN_Q  POST  DONE  FAIL %DONE STATE   DAGNAME                                 
    0     0     0     0     0    26     0 100.0 Success 00/00/level-2/level-2.dag               
    0     0     0     0     0     6     0 100.0 Success *gem-0.dag                              
    0     0     0     0     0    32     0 100.0         TOTALS (32 jobs)                        
Summary: 2 DAGs total (Success:2)

# Final output file is located in
ls -lh runs/osg-gem-1498085153/outputs/merged_GEM.tab 
-rw-r--r-- 1 upendra upendra 68K Jun 21 17:50 runs/osg-gem-1498085153/outputs/merged_GEM.tab 

Monitoring Workflow

Pegasus provides a set of commands to monitor workflow progress. The path to the workflow files as well as commands to monitor the workflow will print to screen upon submitting the workflow. For example:

2016.05.26 23:31:03.859 CDT:   Your workflow has been started and is running in
2016.05.26 23:31:03.869 CDT:     /stash2/user/username/workflows/osg-gem-x/workflow/osg-gem-x
2016.05.26 23:31:03.880 CDT:   *** To monitor the workflow you can run ***
2016.05.26 23:31:03.891 CDT:     pegasus-status -l /stash2/user/username/workflows/osg-gem-x/workflow/osg-gem-x
2016.05.26 23:31:03.901 CDT:   *** To remove your workflow run ***
2016.05.26 23:31:03.912 CDT:     pegasus-remove /stash2/user/username/workflows/osg-gem-x/workflow/osg-gem-x

Output will be transferred at the base of this directory upon completion. For example:

Pre-Workflow User Input

The user must provide indexed reference genome files as well as gene model annotation information prior to submitting the workflow. The user must select a reference prefix($REF_PREFIX) that will be recognized by Pegasus as well as by Hisat2 or Tophat2. In addition, information about splice sites or a reference transcriptome must be provided in order to guide accurate mapping of split input files. Once the user has downloaded a reference genome fasta file and gene annotation in GTF/GFF3 format, the following commands can be used to produce the necessary input files, using GRCh38 as an example $REF_PREFIX for Gencode Release 24 of the human reference genome:

Hisat2:

• Index the reference genome

$ cd reference
$ hisat2-build -f GRCh38.fa GRCh38

• Generate Tab delimited list of splice sites using gene model GTF file as input (Python DefaultDictionary Module necessary)

$ python hisat2_extract_splice_sites.py GRCh38-gencode.v24.annotation.gtf > GRCh38.Splice_Sites.txt

Tophat2:

• Index the reference genome

$ cd reference
$ bowtie2-build -f GRCh38.fa GRCh38

• Generate and Index Reference Transcriptome

$ tophat2 -G GRCh38.gencode.v24.annotation.gff3 --transcriptome-index=transcriptome_data/GRCh38 GRCh38
$ tar czf GRCh38.transcriptome_data.tar.gz transcriptome_data/

STAR:

$ cd reference/star_index
$ bowtie2-build -f GRCh38.fa GRCh38
$ STAR-2.5.2b/bin/Linux_x86_64_static/STAR --runMode genomeGenerate --runThreadN 4 --genomeDir ./ --genomeFastaFiles ../chr21-GRCh38.fa --sjdbGTFfile ../chr21-gencode.v24.annotation.gff3  

• Generate Tab delimited list of splice sites using gene model GTF file as input (Python DefaultDictionary Module necessary)

$ python hisat2_extract_splice_sites.py GRCh38-gencode.v24.annotation.gtf > GRCh38.Splice_Sites.txt