Learn R Programming

GREP2 : GEO RNA-seq Experiments Processing Pipeline

The Gene Expression Omnibus (GEO) is a public repository of gene expression data that hosts more than 6,000 RNA-seq datasets and this number is increasing. Most of these datasets are deposited in raw sequencing format which needs to be downloaded and processed. With an aim to transform all these datasets in an analysis-ready format, we have developed a comprehensive pipeline to simultaneously download and process RNA-seq data sets from GEO. This R-based automated pipeline can process the available RNA-seq data of human, mouse, and rat from GEO. This package is recommended to use in the unix environment as many of the features are not available in windows.


Installation

Before installing GREP2, you need to install the following software packages first:

  1. SRA toolkit
  2. Aspera-connect
  3. FastQC
  4. Salmon
  5. MultiQC

You will also need to install the following R and Bioconductor packages:

install.packages(c("devtools", "XML", "parallel", "utils", "rentrez", "RCurl")
source("https://bioconductor.org/biocLite.R")
biocLite(c("GEOquery", "Biobase", "tximport", "EnsDb.Hsapiens.v86", "EnsDb.Rnorvegicus.v79", "EnsDb.Mmusculus.v79",
    "AnnotationDbi", "org.Hs.eg.db", "org.Mm.eg.db", "org.Rn.eg.db"))

Once you install the above packages, you can now install GREP2 using CRAN:

install.packages("GREP2")

You can also install GREP2 using devtools:

library(devtools)
install_github("uc-bd2k/GREP2")

GREP2 pipeline workflow

To consistently process GEO RNA-seq datasets through a robust and uniform system, we have built GEO RNA-seq evenly processing pipeline (GREP2). To demonstrate the usage of the package, we demonstrate the processing steps with a small dataset from GEO: GSE102170. The whole processing workflow can be summarized in the following steps:

  1. The pipeline starts with a valid GEO series accession ID. Currently

the pipeline works for human, mouse, and rat species. We then retrieve metadata for the GEO series accession using Bioconductor package GEOquery. We also download metadata file from the sequence read archive (SRA) to get corresponding run information.

library(GREP2)
metadata <- get_metadata(geo_series_acc="GSE102170",destdir=tempdir(),
geo_only=FALSE,download_method="auto")
  1. Download corresponding experiment run files from the SRA using ascp

utility of Aspera Connect or regular download. All the downloaded files are stored in the local repository until processed. You can skip this step by downloading fastq files directly.

srr_id <- metadata$metadata_sra$Run
for(i in 1:length(srr_id)){
	get_srr(srr_id=srr_id[i], destdir=tempdir(), ascp=FALSE,
	prefetch_workspace=NULL,ascp_path=NULL)
}
  1. Convert SRA files to fastq format using NCBI SRA toolkit

or download fastq files directly.

library_layout <- metadata$metadata_sra$LibraryLayout
for(i in 1:length(srr_id)){
	get_fastq(srr_id=srr_id[i],library_layout=library_layout[i],
	use_sra_file=FALSE,sra_files_dir=NULL,n_thread=2,
	destdir=tempdir())
}
  1. Run FastQC

on each fastq file to generate quality control (QC) reports.

run_fastqc(destdir=tempdir(),fastq_dir=tempdir(),
n_thread=2)
  1. Remove adapter sequences if necessary using Trimmomatic.
for(i in 1:length(srr_id)){
	trim_fastq(srr_id=srr_id[i],fastq_dir=tempdir(),
	instrument="MiSeq",library_layout=library_layout[i],destdir=tempdir(),n_thread=2)
}
  1. Quantify transcript abundances using Salmon.

Transcript level estimates are then summarized to gene level using Bioconductor package tximport. We obtained gene annotation for Homo sapiens (GRCh38), Mus musculus (GRCm38), and Rattus norvegicus (Rnor_6.0) from Ensemble.

# Before running Salmon, you will have to build index first.
build_index(species="human",kmer=31,ens_release=92,
destdir=tempdir())
# Run Salmon
for(i in 1:length(srr_id)){
	run_salmon(srr_id=srr_id[i],library_layout=library_layout[i],
	index_dir=tempdir(),destdir=tempdir(),
	fastq_dir=tempdir(),use_trimmed_fastq=FALSE,
	other_opts=NULL,n_thread=2)
}
# Run tximport
counts_list <- run_tximport(srr_id=srr_id, species="human",
salmon_dir=paste0(tempdir(),"/salmon"),countsFromAbundance="lengthScaledTPM")
  1. Compile FastQC reports and Salmon log files into a single

interactive HTML report using MultiQC.

run_multiqc(fastqc_dir=tempdir(),salmon_dir=tempdir(),
destdir=tempdir())

You can run the above individual functions for each step or run the whole pipeline using the following process_geo_rnaseq function. All of the above steps are combined into the following single function. We would recommend using this function for processing GEO RNA-seq data.

process_geo_rnaseq (geo_series_acc=geo_series_acc,destdir=tempdir(),
download_method="auto",
ascp=FALSE,prefetch_workspace=NULL,
ascp_path=NULL,use_sra_file=FALSE,trim_fastq=FALSE,
index_dir=tempdir(),species="human",
countsFromAbundance="lengthScaledTPM",n_thread=1)

Copy Link

Version

Install

install.packages('GREP2')

Monthly Downloads

36

Version

1.0.2

License

GPL-3

Issues

Pull Requests

Stars

Forks

Maintainer

Last Published

September 3rd, 2019

Functions in GREP2 (1.0.2)

run_fastqc

QC report for each fastq files using FastQC
run_tximport

Wrapper function to run tximport
trim_fastq

Trim fastq files using Trimmomatic
run_multiqc

Generate combined QC report for Salmon and FastQC
run_salmon

Quantify transcript abundances using Salmon
build_index

Build index for mapping using Salmon
get_metadata

Download metadata from GEO and SRA
get_srr

Download SRA run files
get_fastq

Download fastq files
process_geo_rnaseq

A complete pipeline to process GEO RNA-seq data