An open RNA-Seq data analysis pipeline tutorial with an example of reprocessing data from a recent Zika virus study

Zichen Wang1 and Avi Ma'ayan1*

1Department of Pharmacology and Systems Therapeutics;
BD2K-LINCS Data Coordination and Integration Center;
Mount Sinai Knowledge Management Center for Illuminating the Druggable Genome;
Icahn School of Medicine at Mount Sinai, One Gustave L. Levy Place, Box 1603, New York, NY 10029 USA

*Correspondence: [email protected]


Abstract

RNA-seq analysis is becoming a standard method for global gene expression profiling. However, open and standard pipelines to perform RNA-seq analysis by non-experts remain challenging due to the large size of the raw data files and the hardware requirements for running the alignment step. Here we introduce a reproducible open source RNA-seq pipeline delivered as an IPython notebook and a Docker image. The pipeline uses state-of-the-art tools and can run on various platforms with minimal configuration overhead. The pipeline enables the extraction of knowledge from typical RNA-seq studies by generating interactive principal component analysis (PCA) and hierarchical clustering (HC) plots, performing enrichment analyses against over 90 gene set libraries, and obtaining lists of small molecules that are predicted to either mimic or reverse the observed changes in mRNA expression. We apply the pipeline to a recently published RNA-seq dataset collected from human neuronal progenitors infected with the Zika virus (ZIKV). In addition to confirming the presence of cell cycle genes among the genes that are downregulated by ZIKV, our analysis uncovers significant overlap with upregulated genes that when knocked out in mice induce defects in brain morphology. This result potentially points to the molecular processes associated with the microcephaly phenotype observed in newborns from pregnant mothers infected with the virus. In addition, our analysis predicts small molecules that can either mimic or reverse the expression changes induced by ZIKV. The IPython notebook and Docker image are freely available at: http://nbviewer.jupyter.org/github/maayanlab/Zika-RNAseq-Pipeline/blob/master/Zika.ipynb and https://hub.docker.com/r/maayanlab/zika/

Keywords

Systems biology, bioinformatics pipeline, gene expression analysis, RNA-seq


Introduction

The increase in awareness about the irreproducibility of scientific research requires the development of methods that make experimental and computational protocols easily repeatable and transparent [1]. The advent of interactive notebooks for data analysis pipelines significantly enhances the recording and sharing of data, source code, and figures [2]. In a subset of recent publications, an interactive notebook was published alongside customary manuscripts [3]. Similarly, here we present an interactive IPython notebook (http://nbviewer.jupyter.org/github/maayanlab/Zika-RNAseq-Pipeline/blob/master/Zika.ipynb) that serves as a tutorial for performing a standard RNA-seq pipeline. The IPython notebook pipeline provides scripts (http://dx.doi.org/10.5281/zenodo.56311) that process the raw data into interactive figures and permits other downstream analyses that can enable others to quickly and properly repeat our analysis as well as extract knowledge from their own data. As an example, we applied the pipeline to RNA-seq data from a recent publication where human induced pluripotent stem cells were differentiated to neuronal progenitors and then infected with Zika virus (ZIKV) [4]. The aim of the study was to begin to understand the molecular mechanisms that induce the observed devastating phenotype of newborn-microcephaly from pregnant mothers infected with the virus.

Methods and results

The first publicly available study profiling gene expression changes after ZIKV infection of human cells was deposited into NCBI's Gene Expression Omnibus (GEO) in March 2016. The raw data is available (ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/) from the Sequence Read Archive (SRA) with accession number GSE78711. In this study, gene expression was measured by RNA-seq using two platforms: MiSeq and NextSeq [4] in duplicates. The total number of samples is eight, with four untreated samples and four infected samples. We first downloaded the raw sequencing files from SRA and then converted to FASTQ files. Quality Control (QC) for the RNA-Seq reads was assessed using FastQC [5]. The reports generated by FastQC were in HTML format and can be accessed through hyperlinks from the IPython notebook. The reads in the FASTQ files were aligned to the human genome with Spliced Transcripts Alignment to a Reference (STAR) [6]. STAR is a leading aligner that accomplishes the alignment step faster and more accurately than other current alternatives [6]. We next applied featureCounts [7] to assign reads to genes, and then applied the edgeR Bioconductor package [8] to compute counts per million (CPM) and reads per kilobase million (RPKM). The next steps are performed in Python within the IPython notebook. We first filtered out genes that are not expressed or lowly expressed. Subsequently, we performed principal component analysis (PCA) (Fig. 1). The PCA plots show that the samples cluster by infected vs. control cells, but also by platform. Next, we visualized the 800 genes with the largest variance using an interactive hierarchical clustering (HC) plot (Fig. 2). This analysis separates the groups of genes that are differentially expressed by infected vs. control from those that are differential by platform. The visualization of the clusters is implemented with an interactive external web-based data visualization tool called clustergrammer (http://amp.pharm.mssm.edu/clustergrammer/). Clustergrammer provides interactive searching, sorting and zoom capabilities.

The following step is to identify the differentially expressed genes (DEG) between the two conditions. This is achieved with a unique method we developed called the Characteristic Direction (CD) [9]. The CD method is a multivariate method that we have previously demonstrated to outperform other leading methods that compute differential expression between two conditions [9]. Once we have ranked the lists of DEG, we submit these for signature analysis using two tools: Enrichr [10] and L1000CDS2 [11]. Enrichr queries the up and down gene sets against over 180,000 annotated gene sets belonging to 90 gene set libraries covering pathway databases, ontologies, disease databases, and more [10]. The results from this enrichment analysis confirm that the downregulated genes after ZIKV infection are enriched for genes involved in cell cycle-related processes (Fig. 3a). These genes are enriched for targets of the transcription factors E2F4 and FOXM1 (Fig. 3b). Both transcription factors are known to regulate cell proliferation and play central role in many cancers. The downregulation of cell cycle genes was already reported in the original publication; nevertheless, we obtained more interesting results for the enriched terms that appeared most significant for the upregulated genes. Particularly, the top two terms from the mouse genome informatics (MGI) Mammalian Phenotype Level 4 library are abnormal nervous system (MP0003861) and abnormal brain morphology (MP0002152) (Table S1). This library associates gene knockouts in mice with mammalian phenotypes. These enriched terms enlist a short set of genes that potentially link ZIKV infection with the concerning observed microcephaly phenotype. Finally, to identify small molecules that can potentially either reverse or mimic ZIKV-induced gene expression changes, we query the ZIKV-induced signatures against the LINCS L1000 data. For this, we utilize L1000CDS2 [11], a search engine that prioritize small molecules given a gene expression signature as input. L1000CDS2 contains 30,000 significant signatures that were processed from the LINCS L1000 data with the CD method. The results suggest small molecules that could be tested in follow-up studies in human cells for potential efficacy against ZIKV (Table S2).

To ensure the reproducibility of the computational environment used for the whole RNA-Seq pipeline, we packaged all the software components used in this tutorial, including the command line tools, R packages, and Python packages into a Docker image. This Docker image is made publically available at https://hub.docker.com/r/maayanlab/zika/. The Docker image was created based on the specifications outlined on the official IPython’s Scipy Stack image (https://hub.docker.com/r/IPython/scipystack/). The additional command line tools, R scripts, and Python packages together with their dependencies were compiled and installed into the Docker image. The RNA-Seq pipeline Docker image was deployed onto our Mesos cluster, which allows users to run the IPython notebook interactively. The Docker image can also be downloaded and executed on local computers and servers, or deployed in the cloud if users have access to cloud provider services with a Docker Toolbox installed (https://www.docker.com/products/docker-toolbox). We also provide detailed instructions on how to download and execute the Docker image (https://hub.docker.com/r/maayanlab/zika/).

The ‘Dockerization’ of the RNA-Seq pipeline facilitates reproducibility of the pipeline at the software level because the Docker image ensures that all versions of the software components are consistent and static. Dockerization also helps users to handle the complex installation of many dependencies required for the computational pipeline. Moreover, the Docker image can be executed on a single computer, clusters/servers and on the cloud. The only limitation of having a Docker image is that it prevents users from adding or altering the various steps which require additional software components and packages. However, advanced users can build their own Docker images based on our initial image to customize it for their needs.

In [1]:
import os
import numpy as np
import pandas as pd

Below we assign some global variables that will be used across the rest of the notebook.
Please change these variables accordingly if you intend to use this for other studies.

In [2]:
# The URL for the SRA study (project), usually in a SRPxxxxx folder including several SRRxxxxx folders (samples)
os.environ['FTP_URL'] = 'ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/'
# The working directory to store all the sequencing data, will be created if not exists
os.environ['WORKDIR'] = 'data/Zika/'
# The directory for the reference genome
os.environ['GENOMEDIR'] = 'genomes/Homo_sapiens/UCSC/hg19'
In [3]:
# Download the SRA files
# This bash script is commented out because we don't want to download the files every time.
!mkdir -p $WORKDIR
!wget -r $FTP_URL \
    --no-parent -nH --cut-dirs=8 \
    -P $WORKDIR
--2016-06-23 11:53:28--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/
           => 'data/Zika/.listing'
Resolving ftp-trace.ncbi.nlm.nih.gov... 130.14.250.13, 2607:f220:41e:250::7
Connecting to ftp-trace.ncbi.nlm.nih.gov|130.14.250.13|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895 ... done.
==> PASV ... done.    ==> LIST ... done.

.listing                [ <=>                ]     600  --.-KB/s    in 0.003s  

2016-06-23 11:53:28 (209 KB/s) - 'data/Zika/.listing' saved [600]

Removed 'data/Zika/.listing'.
--2016-06-23 11:53:28--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191542/
           => 'data/Zika/SRR3191542/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191542 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3191542/.listing     [ <=>                ]      73  --.-KB/s    in 0.001s  

2016-06-23 11:53:28 (84.3 KB/s) - 'data/Zika/SRR3191542/.listing' saved [73]

Removed 'data/Zika/SRR3191542/.listing'.
--2016-06-23 11:53:28--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191542/SRR3191542.sra
           => 'data/Zika/SRR3191542/SRR3191542.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3191542.sra ... done.
Length: 486949093 (464M)

SRR3191542/SRR31915 100%[===================>] 464.39M  84.8MB/s    in 6.0s    

2016-06-23 11:53:34 (77.8 MB/s) - 'data/Zika/SRR3191542/SRR3191542.sra' saved [486949093]

--2016-06-23 11:53:34--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191543/
           => 'data/Zika/SRR3191543/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191543 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3191543/.listing     [ <=>                ]      73  --.-KB/s    in 0.001s  

2016-06-23 11:53:34 (63.3 KB/s) - 'data/Zika/SRR3191543/.listing' saved [73]

Removed 'data/Zika/SRR3191543/.listing'.
--2016-06-23 11:53:34--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191543/SRR3191543.sra
           => 'data/Zika/SRR3191543/SRR3191543.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3191543.sra ... done.
Length: 465054386 (444M)

SRR3191543/SRR31915 100%[===================>] 443.51M  93.4MB/s    in 5.5s    

2016-06-23 11:53:40 (80.4 MB/s) - 'data/Zika/SRR3191543/SRR3191543.sra' saved [465054386]

--2016-06-23 11:53:40--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191544/
           => 'data/Zika/SRR3191544/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191544 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3191544/.listing     [ <=>                ]      73  --.-KB/s    in 0.003s  

2016-06-23 11:53:40 (25.4 KB/s) - 'data/Zika/SRR3191544/.listing' saved [73]

Removed 'data/Zika/SRR3191544/.listing'.
--2016-06-23 11:53:40--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191544/SRR3191544.sra
           => 'data/Zika/SRR3191544/SRR3191544.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3191544.sra ... done.
Length: 450651106 (430M)

SRR3191544/SRR31915 100%[===================>] 429.77M   101MB/s    in 5.0s    

2016-06-23 11:53:45 (86.7 MB/s) - 'data/Zika/SRR3191544/SRR3191544.sra' saved [450651106]

--2016-06-23 11:53:45--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191545/
           => 'data/Zika/SRR3191545/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191545 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3191545/.listing     [ <=>                ]      73  --.-KB/s    in 0.002s  

2016-06-23 11:53:45 (39.2 KB/s) - 'data/Zika/SRR3191545/.listing' saved [73]

Removed 'data/Zika/SRR3191545/.listing'.
--2016-06-23 11:53:45--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3191545/SRR3191545.sra
           => 'data/Zika/SRR3191545/SRR3191545.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3191545.sra ... done.
Length: 467255848 (446M)

SRR3191545/SRR31915 100%[===================>] 445.61M   103MB/s    in 4.9s    

2016-06-23 11:53:50 (91.8 MB/s) - 'data/Zika/SRR3191545/SRR3191545.sra' saved [467255848]

--2016-06-23 11:53:50--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194428/
           => 'data/Zika/SRR3194428/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194428 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3194428/.listing     [ <=>                ]      74  --.-KB/s    in 0.002s  

2016-06-23 11:53:50 (33.3 KB/s) - 'data/Zika/SRR3194428/.listing' saved [74]

Removed 'data/Zika/SRR3194428/.listing'.
--2016-06-23 11:53:50--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194428/SRR3194428.sra
           => 'data/Zika/SRR3194428/SRR3194428.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3194428.sra ... done.
Length: 2499769558 (2.3G)

SRR3194428/SRR31944 100%[===================>]   2.33G  89.4MB/s    in 27s     

2016-06-23 11:54:17 (88.0 MB/s) - 'data/Zika/SRR3194428/SRR3194428.sra' saved [2499769558]

--2016-06-23 11:54:17--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194429/
           => 'data/Zika/SRR3194429/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194429 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3194429/.listing     [ <=>                ]      74  --.-KB/s    in 0.001s  

2016-06-23 11:54:17 (131 KB/s) - 'data/Zika/SRR3194429/.listing' saved [74]

Removed 'data/Zika/SRR3194429/.listing'.
--2016-06-23 11:54:17--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194429/SRR3194429.sra
           => 'data/Zika/SRR3194429/SRR3194429.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3194429.sra ... done.
Length: 3242990685 (3.0G)

SRR3194429/SRR31944 100%[===================>]   3.02G  95.5MB/s    in 35s     

2016-06-23 11:54:52 (89.2 MB/s) - 'data/Zika/SRR3194429/SRR3194429.sra' saved [3242990685]

--2016-06-23 11:54:52--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194430/
           => 'data/Zika/SRR3194430/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194430 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3194430/.listing     [ <=>                ]      74  --.-KB/s    in 0.001s  

2016-06-23 11:54:52 (132 KB/s) - 'data/Zika/SRR3194430/.listing' saved [74]

Removed 'data/Zika/SRR3194430/.listing'.
--2016-06-23 11:54:52--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194430/SRR3194430.sra
           => 'data/Zika/SRR3194430/SRR3194430.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3194430.sra ... done.
Length: 2616536610 (2.4G)

SRR3194430/SRR31944 100%[===================>]   2.44G  81.9MB/s    in 31s     

2016-06-23 11:55:23 (81.0 MB/s) - 'data/Zika/SRR3194430/SRR3194430.sra' saved [2616536610]

--2016-06-23 11:55:23--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194431/
           => 'data/Zika/SRR3194431/.listing'
==> CWD (1) /sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194431 ... done.
==> PASV ... done.    ==> LIST ... done.

SRR3194431/.listing     [ <=>                ]      74  --.-KB/s    in 0.001s  

2016-06-23 11:55:23 (132 KB/s) - 'data/Zika/SRR3194431/.listing' saved [74]

Removed 'data/Zika/SRR3194431/.listing'.
--2016-06-23 11:55:23--  ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP070/SRP070895/SRR3194431/SRR3194431.sra
           => 'data/Zika/SRR3194431/SRR3194431.sra'
==> CWD not required.
==> PASV ... done.    ==> RETR SRR3194431.sra ... done.
Length: 2287601650 (2.1G)

SRR3194431/SRR31944 100%[===================>]   2.13G  94.4MB/s    in 24s     

2016-06-23 11:55:47 (89.5 MB/s) - 'data/Zika/SRR3194431/SRR3194431.sra' saved [2287601650]

FINISHED --2016-06-23 11:55:47--
Total wall clock time: 2m 20s
Downloaded: 8 files, 12G in 2m 18s (86.4 MB/s)
In [4]:
## Examine the downloaded SRA files downloaded
!ls -lh $WORKDIR/SRR*/*.sra
-rw-r--r--  1 zichen  staff   464M Feb 26 21:21 data/Zika//SRR3191542/SRR3191542.sra
-rw-r--r--  1 zichen  staff   444M Feb 26 21:22 data/Zika//SRR3191543/SRR3191543.sra
-rw-r--r--  1 zichen  staff   430M Feb 26 21:21 data/Zika//SRR3191544/SRR3191544.sra
-rw-r--r--  1 zichen  staff   446M Feb 26 21:21 data/Zika//SRR3191545/SRR3191545.sra
-rw-r--r--  1 zichen  staff   2.3G Feb 29 17:19 data/Zika//SRR3194428/SRR3194428.sra
-rw-r--r--  1 zichen  staff   3.0G Feb 29 17:20 data/Zika//SRR3194429/SRR3194429.sra
-rw-r--r--  1 zichen  staff   2.4G Feb 29 17:21 data/Zika//SRR3194430/SRR3194430.sra
-rw-r--r--  1 zichen  staff   2.1G Feb 29 17:09 data/Zika//SRR3194431/SRR3194431.sra

The downloaded SRA files are next processed by following these steps:

  1. fastq-dump in the SRA-toolkit to generate .fastq files
  2. FastQC[3](#ref3) to perform Quality Controls and generate QC report for the input RNA-seq data
  3. STAR[3](#ref3) for the read alignment
  4. featureCounts[4](#ref4) for assigning reads to genes
  5. edgeR Bioconductor package[5](#ref5) were used to compute CPM and RPKM

Steps 1-4 are processed by this bash script analyze_sra.sh. This bash script can take command line arguments specifying the location of the reference genome and working directory for the SRA files.

In [5]:
!bash analyze_sra.sh -h
Usage: ./analyze_sra.sh -g <GENOME> -w <WORKDIR>

We run the bash script by specifying the working directory and the genome directory and pipe the log into a analyze_sra.log file

In [6]:
!bash analyze_sra.sh -w $WORKDIR -g $GENOMEDIR | tee analyze_sra.log
Dumping .sra to .fastq
SINGLE/SRR3194428/SRR3194428.sra
Read 72983243 spots for SINGLE/SRR3194428/SRR3194428.sra
Written 72983243 spots for SINGLE/SRR3194428/SRR3194428.sra
SINGLE/SRR3194429/SRR3194429.sra
Read 94729809 spots for SINGLE/SRR3194429/SRR3194429.sra
Written 94729809 spots for SINGLE/SRR3194429/SRR3194429.sra
SINGLE/SRR3194430/SRR3194430.sra
Read 76299868 spots for SINGLE/SRR3194430/SRR3194430.sra
Written 76299868 spots for SINGLE/SRR3194430/SRR3194430.sra
SINGLE/SRR3194431/SRR3194431.sra
Read 66528035 spots for SINGLE/SRR3194431/SRR3194431.sra
Written 66528035 spots for SINGLE/SRR3194431/SRR3194431.sra
PAIRED/SRR3191542/SRR3191542.sra
Read 7927777 spots for PAIRED/SRR3191542/SRR3191542.sra
Written 7927777 spots for PAIRED/SRR3191542/SRR3191542.sra
PAIRED/SRR3191543/SRR3191543.sra
Read 7391076 spots for PAIRED/SRR3191543/SRR3191543.sra
Written 7391076 spots for PAIRED/SRR3191543/SRR3191543.sra
PAIRED/SRR3191544/SRR3191544.sra
Read 7361527 spots for PAIRED/SRR3191544/SRR3191544.sra
Written 7361527 spots for PAIRED/SRR3191544/SRR3191544.sra
PAIRED/SRR3191545/SRR3191545.sra
Read 7621347 spots for PAIRED/SRR3191545/SRR3191545.sra
Written 7621347 spots for PAIRED/SRR3191545/SRR3191545.sra
Started to align reads to the genome and assemble transcriptome
SRR3194428
Started analysis of SRR3194428.fastq
Approx 5% complete for SRR3194428.fastq
Approx 10% complete for SRR3194428.fastq
Approx 15% complete for SRR3194428.fastq
Approx 20% complete for SRR3194428.fastq
Approx 25% complete for SRR3194428.fastq
Approx 30% complete for SRR3194428.fastq
Approx 35% complete for SRR3194428.fastq
Approx 40% complete for SRR3194428.fastq
Approx 45% complete for SRR3194428.fastq
Approx 50% complete for SRR3194428.fastq
Approx 55% complete for SRR3194428.fastq
Approx 60% complete for SRR3194428.fastq
Approx 65% complete for SRR3194428.fastq
Approx 70% complete for SRR3194428.fastq
Approx 75% complete for SRR3194428.fastq
Approx 80% complete for SRR3194428.fastq
Approx 85% complete for SRR3194428.fastq
Approx 90% complete for SRR3194428.fastq
Approx 95% complete for SRR3194428.fastq
Analysis complete for SRR3194428.fastq
Jun 15 17:29:04 ..... Started STAR run
Jun 15 17:29:04 ..... Loading genome
Jun 15 17:30:08 ..... Processing annotations GTF
Jun 15 17:30:12 ..... Inserting junctions into the genome indices
Jun 15 17:31:38 ..... Started mapping
Jun 15 17:37:54 ..... Started sorting BAM
Jun 15 17:40:11 ..... Finished successfully

        ==========     _____ _    _ ____  _____  ______          _____  
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
	v1.4.6-p2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S ../star_output/SRR3194428Aligned.sortedByC ... ||
||                                                                            ||
||             Output file : ../featureCount_output/SRR3194428.count.txt      ||
||             Annotations : /notebook/genomes/Homo_sapiens/UCSC/hg19/Ann ... ||
||                                                                            ||
||                 Threads : 8                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file /notebook/genomes/Homo_sapiens/UCSC/hg19/Annotati ... ||

Warning: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id' 
The attributes included in your GTF annotation are 'gene_id ""; gene_name ""; transcript_id "NR_001526_1"; tss_id "TSS25560";' 

||    Features : 482184                                                       ||
||    Meta-features : 25370                                                   ||
||    Chromosomes : 51                                                        ||
||                                                                            ||
|| Process BAM file ../star_output/SRR3194428Aligned.sortedByCoord.out.ba ... ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 78482481                                                  ||
||    Successfully assigned reads : 48156708 (61.4%)                          ||
||    Running time : 1.06 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

SRR3194429
Started analysis of SRR3194429.fastq
Approx 5% complete for SRR3194429.fastq
Approx 10% complete for SRR3194429.fastq

Approx 20% complete for SRR3194429.fastq
Approx 25% complete for SRR3194429.fastq
Approx 30% complete for SRR3194429.fastq
Approx 35% complete for SRR3194429.fastq
Approx 40% complete for SRR3194429.fastq
Approx 45% complete for SRR3194429.fastq
Approx 50% complete for SRR3194429.fastq
Approx 55% complete for SRR3194429.fastq
Approx 60% complete for SRR3194429.fastq
Approx 65% complete for SRR3194429.fastq
Approx 70% complete for SRR3194429.fastq
Approx 75% complete for SRR3194429.fastq
Approx 80% complete for SRR3194429.fastq
Approx 85% complete for SRR3194429.fastq
Approx 90% complete for SRR3194429.fastq
Approx 95% complete for SRR3194429.fastq
Analysis complete for SRR3194429.fastq
Jun 15 17:48:09 ..... Started STAR run
Jun 15 17:48:10 ..... Loading genome
Jun 15 17:49:13 ..... Processing annotations GTF
Jun 15 17:49:17 ..... Inserting junctions into the genome indices
Jun 15 17:50:43 ..... Started mapping
Jun 15 17:58:51 ..... Started sorting BAM
Jun 15 18:01:58 ..... Finished successfully

        ==========     _____ _    _ ____  _____  ______          _____  
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
	v1.4.6-p2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S ../star_output/SRR3194429Aligned.sortedByC ... ||
||                                                                            ||
||             Output file : ../featureCount_output/SRR3194429.count.txt      ||
||             Annotations : /notebook/genomes/Homo_sapiens/UCSC/hg19/Ann ... ||
||                                                                            ||
||                 Threads : 8                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file /notebook/genomes/Homo_sapiens/UCSC/hg19/Annotati ... ||

Warning: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id' 
The attributes included in your GTF annotation are 'gene_id ""; gene_name ""; transcript_id "NR_001526_1"; tss_id "TSS25560";' 

||    Features : 482184                                                       ||
||    Meta-features : 25370                                                   ||
||    Chromosomes : 51                                                        ||
||                                                                            ||
|| Process BAM file ../star_output/SRR3194429Aligned.sortedByCoord.out.ba ... ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 101843560                                                 ||
||    Successfully assigned reads : 63286307 (62.1%)                          ||
||    Running time : 1.38 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

SRR3194430
Started analysis of SRR3194430.fastq
Approx 5% complete for SRR3194430.fastq
Approx 10% complete for SRR3194430.fastq
Approx 15% complete for SRR3194430.fastq
Approx 20% complete for SRR3194430.fastq
Approx 25% complete for SRR3194430.fastq
Approx 30% complete for SRR3194430.fastq
Approx 35% complete for SRR3194430.fastq
Approx 40% complete for SRR3194430.fastq
Approx 45% complete for SRR3194430.fastq
Approx 50% complete for SRR3194430.fastq
Approx 55% complete for SRR3194430.fastq
Approx 60% complete for SRR3194430.fastq
Approx 65% complete for SRR3194430.fastq
Approx 70% complete for SRR3194430.fastq
Approx 75% complete for SRR3194430.fastq
Approx 80% complete for SRR3194430.fastq
Approx 85% complete for SRR3194430.fastq
Approx 90% complete for SRR3194430.fastq
Approx 95% complete for SRR3194430.fastq
Analysis complete for SRR3194430.fastq
Jun 15 18:08:52 ..... Started STAR run
Jun 15 18:08:52 ..... Loading genome
Jun 15 18:09:54 ..... Processing annotations GTF
Jun 15 18:09:59 ..... Inserting junctions into the genome indices
Jun 15 18:11:27 ..... Started mapping
Jun 15 18:18:06 ..... Started sorting BAM
Jun 15 18:20:37 ..... Finished successfully

        ==========     _____ _    _ ____  _____  ______          _____  
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
	v1.4.6-p2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S ../star_output/SRR3194430Aligned.sortedByC ... ||
||                                                                            ||
||             Output file : ../featureCount_output/SRR3194430.count.txt      ||
||             Annotations : /notebook/genomes/Homo_sapiens/UCSC/hg19/Ann ... ||
||                                                                            ||
||                 Threads : 8                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file /notebook/genomes/Homo_sapiens/UCSC/hg19/Annotati ... ||

Warning: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id' 
The attributes included in your GTF annotation are 'gene_id ""; gene_name ""; transcript_id "NR_001526_1"; tss_id "TSS25560";' 

||    Features : 482184                                                       ||
||    Meta-features : 25370                                                   ||
||    Chromosomes : 51                                                        ||
||                                                                            ||
|| Process BAM file ../star_output/SRR3194430Aligned.sortedByCoord.out.ba ... ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 80847598                                                  ||
||    Successfully assigned reads : 49763125 (61.6%)                          ||
||    Running time : 1.06 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

SRR3194431
Started analysis of SRR3194431.fastq
Approx 5% complete for SRR3194431.fastq
Approx 10% complete for SRR3194431.fastq
Approx 15% complete for SRR3194431.fastq
Approx 20% complete for SRR3194431.fastq
Approx 25% complete for SRR3194431.fastq
Approx 30% complete for SRR3194431.fastq
Approx 35% complete for SRR3194431.fastq
Approx 40% complete for SRR3194431.fastq
Approx 45% complete for SRR3194431.fastq
Approx 50% complete for SRR3194431.fastq
Approx 55% complete for SRR3194431.fastq
Approx 60% complete for SRR3194431.fastq
Approx 65% complete for SRR3194431.fastq
Approx 70% complete for SRR3194431.fastq
Approx 75% complete for SRR3194431.fastq
Approx 80% complete for SRR3194431.fastq
Approx 85% complete for SRR3194431.fastq
Approx 90% complete for SRR3194431.fastq
Approx 95% complete for SRR3194431.fastq
Analysis complete for SRR3194431.fastq
Jun 15 18:26:29 ..... Started STAR run
Jun 15 18:26:29 ..... Loading genome
Jun 15 18:27:31 ..... Processing annotations GTF
Jun 15 18:27:35 ..... Inserting junctions into the genome indices
Jun 15 18:28:59 ..... Started mapping
Jun 15 18:34:53 ..... Started sorting BAM
Jun 15 18:37:00 ..... Finished successfully

        ==========     _____ _    _ ____  _____  ______          _____  
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
	v1.4.6-p2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S ../star_output/SRR3194431Aligned.sortedByC ... ||
||                                                                            ||
||             Output file : ../featureCount_output/SRR3194431.count.txt      ||
||             Annotations : /notebook/genomes/Homo_sapiens/UCSC/hg19/Ann ... ||
||                                                                            ||
||                 Threads : 8                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file /notebook/genomes/Homo_sapiens/UCSC/hg19/Annotati ... ||

Warning: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id' 
The attributes included in your GTF annotation are 'gene_id ""; gene_name ""; transcript_id "NR_001526_1"; tss_id "TSS25560";' 

||    Features : 482184                                                       ||
||    Meta-features : 25370                                                   ||
||    Chromosomes : 51                                                        ||
||                                                                            ||
|| Process BAM file ../star_output/SRR3194431Aligned.sortedByCoord.out.ba ... ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 70336033                                                  ||
||    Successfully assigned reads : 43391323 (61.7%)                          ||
||    Running time : 0.96 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

Step 5 is done with this R script normalize.R

In [7]:
!Rscript normalize.R $WORKDIR
Setting WORKDIR to: data/Zika/ 
Loading required package: limma
Loading required package: methods
Finished reading featureCount into data.frame with shape:  25370 9
RPKM matrix written to "repRpkmMatrix_featureCounts.csv"
CPM matrix file written to "repCpmMatrix_featureCounts.csv"
In [8]:
## We can examine the QC reports from the FastQC program to evaluate the quality of the data
from IPython.display import FileLinks
FileLinks(os.path.join(os.environ['WORKDIR'], 'fastQC_output'), included_suffixes=['.html'])
In [9]:
## Check the alignment stats
## This will output the first 10 lines of all summary files from the featureCounts folder
!head $WORKDIR/featureCount_output/*.summary 
==> data/Zika//featureCount_output/SRR3191542.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3191542Aligned.sortedByCoord.out.bam
Assigned	10859768
Unassigned_Ambiguity	277433
Unassigned_MultiMapping	1315262
Unassigned_NoFeatures	3658745
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3191543.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3191543Aligned.sortedByCoord.out.bam
Assigned	10167142
Unassigned_Ambiguity	264366
Unassigned_MultiMapping	1229314
Unassigned_NoFeatures	3217017
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3191544.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3191544Aligned.sortedByCoord.out.bam
Assigned	10024864
Unassigned_Ambiguity	243520
Unassigned_MultiMapping	1148955
Unassigned_NoFeatures	3421617
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3191545.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3191545Aligned.sortedByCoord.out.bam
Assigned	10340131
Unassigned_Ambiguity	248114
Unassigned_MultiMapping	1177805
Unassigned_NoFeatures	3531209
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3194428.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3194428Aligned.sortedByCoord.out.bam
Assigned	48156708
Unassigned_Ambiguity	1115286
Unassigned_MultiMapping	12147267
Unassigned_NoFeatures	17063220
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3194429.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3194429Aligned.sortedByCoord.out.bam
Assigned	63286307
Unassigned_Ambiguity	1498833
Unassigned_MultiMapping	15919347
Unassigned_NoFeatures	21139073
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3194430.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3194430Aligned.sortedByCoord.out.bam
Assigned	49763125
Unassigned_Ambiguity	1096905
Unassigned_MultiMapping	11709443
Unassigned_NoFeatures	18278125
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

==> data/Zika//featureCount_output/SRR3194431.count.txt.summary <==
Status	/home/maayanlab/Zika/star_output/SRR3194431Aligned.sortedByCoord.out.bam
Assigned	43391323
Unassigned_Ambiguity	952301
Unassigned_MultiMapping	10127586
Unassigned_NoFeatures	15864823
Unassigned_Unmapped	0
Unassigned_MappingQuality	0
Unassigned_FragementLength	0
Unassigned_Chimera	0
Unassigned_Secondary	0

After you completed successfully the above steps, you can start to analyze the processed expression matrix of gene expression in Python

In [10]:
## Load the expression matrix
expr_df = pd.read_csv(os.path.join(os.environ['WORKDIR'], 'repCpmMatrix_featureCounts.csv'))
expr_df = expr_df.set_index(expr_df.columns[0])
expr_df.head()
Out[10]:
SRR3191542 SRR3191543 SRR3191544 SRR3191545 SRR3194428 SRR3194429 SRR3194430 SRR3194431
Unnamed: 0
DDX11L1 0.000000 0.000000 0.199504 0.00000 0.000000 0.000000 0.020095 0.000000
WASH7P 16.759106 18.097515 17.057588 23.59738 2.491865 2.907422 2.190377 2.396793
MIR6859-2 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000
MIR6859-1 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000
FAM138A 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000
In [11]:
print expr_df.shape
(25370, 8)
In [12]:
## Filter out non-expressed genes
expr_df = expr_df.loc[expr_df.sum(axis=1) > 0, :]
print expr_df.shape

## Filter out lowly expressed genes
mask_low_vals = (expr_df > 0.3).sum(axis=1) > 2
expr_df = expr_df.loc[mask_low_vals, :]

print expr_df.shape
(21983, 8)
(16242, 8)

Obtain more metadata about the samples by clicking on the RunInfo Table button on the SRP page available in this URL http://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP070895. Clicking on this button downloads a spreadsheet with additional metadata. Next read this file and extract the relavant variables from it.

In [13]:
meta_df = pd.read_csv(os.path.join(os.environ['WORKDIR'], 'SraRunTable.txt'), sep='\t').set_index('Run_s')
print meta_df.shape
# re-order the index to make it the same with expr_df
meta_df = meta_df.ix[expr_df.columns]
meta_df
(8, 27)
Out[13]:
BioSample_s Experiment_s LibraryLayout_s LoadDate_s MBases_l MBytes_l SRA_Sample_s Sample_Name_s infection_status_s Assay_Type_s ... Library_Name_s Organism_s Platform_s ReleaseDate_s SRA_Study_s cell_type_s g1k_analysis_group_s g1k_pop_code_s source_s source_name_s
SRR3191542 SAMN04517925 SRX1602854 PAIRED 2016-02-26 1141 464 SRS1312741 GSM2073121 mock infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3191543 SAMN04517926 SRX1602855 PAIRED 2016-02-26 1063 443 SRS1312740 GSM2073122 mock infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3191544 SAMN04517927 SRX1602856 PAIRED 2016-02-26 1059 429 SRS1312738 GSM2073123 Zika infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3191545 SAMN04517928 SRX1602857 PAIRED 2016-02-26 1095 445 SRS1312739 GSM2073124 Zika infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3194428 SAMN04521221 SRX1605077 SINGLE 2016-02-29 5248 2383 SRS1314803 GSM2075585 mock infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3194429 SAMN04521222 SRX1605078 SINGLE 2016-02-29 6806 3092 SRS1314802 GSM2075586 mock infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3194430 SAMN04521223 SRX1605079 SINGLE 2016-02-29 5481 2495 SRS1314801 GSM2075587 Zika infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells
SRR3194431 SAMN04521224 SRX1605080 SINGLE 2016-02-29 4776 2181 SRS1314800 GSM2075588 Zika infected RNA-Seq ... <not provided> Homo sapiens ILLUMINA 2016-03-04 SRP070895 Neural Progenitor cells <not provided> <not provided> <not provided> human Neural Progenitor cells

8 rows × 27 columns

Now we have everything setup, the first thing to do is to generate PCA plots to observe whether the samples cluster as expected: controls with controls, and treatments with treatments.

In [14]:
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['pdf.fonttype'] = 42 ## Output Type 3 (Type3) or Type 42 (TrueType)
rcParams['font.sans-serif'] = 'Arial'
# ignore FutureWarning that may pop up when plotting
import warnings
warnings.filterwarnings("ignore")
import urllib3
urllib3.disable_warnings()
In [15]:
from IPython.display import HTML, display
# to display hyperlink as <a> tag in output cells
def display_link(url):
    raw_html = '<a href="%s" target="_blank">%s</a>' % (url, url)
    return display(HTML(raw_html))

You can obtain the script RNAseq from this repo.

In [16]:
import RNAseq
In [17]:
# plot PCA
%matplotlib inline
RNAseq.PCA_plot(expr_df.values, meta_df['infection_status_s'], 
         standardize=2, log=True, 
         show_text=False, sep=' ', legend_loc='upper right')

The PCA plot below is the same as above, except that we color the samples by platform (two Illumina sequencing machines used: Illumina MiSeq for paired-end and Illumina NextSeq 500 for single-end).

In [18]:
d_layout_platform = {'PAIRED': 'MiSeq', 'SINGLE': 'NextSeq 500'}
meta_df['platform'] = [d_layout_platform[l] for l in meta_df['LibraryLayout_s']]
RNAseq.PCA_plot(expr_df.values, meta_df['platform'], 
         standardize=2, log=True, 
         show_text=False, sep=' ', legend_loc='upper right')

We can also plot a 3D interactive PCA plot using plotly, which has a nice integration with Jupyter notebooks.

In [19]:
# Compute the coordinates of samples in the PCA space
variance_explained, pca_transformed = RNAseq.perform_PCA(expr_df.values, standardize=2, log=True)
# Bind x, y, z coordinates to meta_df
meta_df['x'] = pca_transformed[:,0]
meta_df['y'] = pca_transformed[:,1]
meta_df['z'] = pca_transformed[:,2]
In [20]:
import plotly
plotly.offline.init_notebook_mode() # To embed plots in the output cell of the notebook
import plotly.graph_objs as go

conditions = meta_df['infection_status_s'].unique().tolist()
platforms = meta_df['platform'].unique().tolist()
SYMBOLS = ['circle', 'square']
COLORS = RNAseq.COLORS10

data = [] # To collect all Scatter3d instances
for (condition, platform), meta_df_sub in meta_df.groupby(['infection_status_s', 'platform']):
    # Iteratate through samples grouped by condition and platform
    display_name = '%s, %s' % (condition, platform)
    # Initiate a Scatter3d instance for each group of samples specifying their coordinates
    # and displaying attributes including color, shape, size and etc.
    trace = go.Scatter3d(
        x=meta_df_sub['x'],
        y=meta_df_sub['y'],
        z=meta_df_sub['z'],
        text=meta_df_sub.index,
        mode='markers',
        marker=dict(
            size=10,
            color=COLORS[conditions.index(condition)], # Color by infection status
            symbol=SYMBOLS[platforms.index(platform)], # Shaped by sequencing platforms
            opacity=.8,
        ),
        name=display_name,
    )
    
    data.append(trace)

# Configs for layout and axes
layout=dict(height=1000, width=1000, 
            title='3D PCA plot for samples in Zika study',
            scene=dict(
                xaxis=dict(title='PC1 (%.2f%% variance)' % variance_explained[0]),
                yaxis=dict(title='PC2 (%.2f%% variance)' % variance_explained[1]),
                zaxis=dict(title='PC3 (%.2f%% variance)' % variance_explained[2])
               )
)
fig=dict(data=data, layout=layout)
plotly.offline.iplot(fig)
/Library/Python/2.7/site-packages/requests/packages/urllib3/util/ssl_.py:315: SNIMissingWarning: An HTTPS request has been made, but the SNI (Subject Name Indication) extension to TLS is not available on this platform. This may cause the server to present an incorrect TLS certificate, which can cause validation failures. For more information, see https://urllib3.readthedocs.org/en/latest/security.html#snimissingwarning.
  SNIMissingWarning