Next Generation Sequencing RNAseq Data Analysis (Part 1)

Preface

Here I am going to write a series of posts dedicated to high throughput RNA sequencing data analysis. As an essential technology, RNA sequencing is applied to many fields of biology and life science. Knowing how to process RNA seq data is important for researchers in these fields. This series of articles will explain the steps involved in the process.

In our first article, we will obtain real-world Next Generation Sequencing (NGS) data, learn about the data format required for further analysis, and annotate the data accordingly.

Where to get a data

For this project we will get data from NCBI Gene Expression Omnibus Database (GEO) which is a public repository of functional genomics data. This database contains thousands of datasets related to analysis of gene expression. We will download one of them (Series GSE193895). Results of this research were published here¹.

In series of my articles I’m not going to replicate the original research but describe the ways to analyse NGS RNA seq data.

Downloading data

We need to download the archive with raw counts for each sample GSE193895_RAW.tar. It contains 35 separate files named after sample id followed by treatment conditions. Extract them tar -xf GSE193895_RAW.tar. The files are also archived so use gzip -dv * to unarchive all of them.

Data wrangling

Load all needed packages:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

# run explicit error when packages contain functions with same names
library(conflicted)

# for easier data manipulations
library(tidyverse)

# for pretty table printing
library(knitr)

# for splitting strings
library(stringr)

# for annotation
library(AnnotationHub)

library(ensembldb)

Get list of all files with raw counts for each sample:

1

count_files <- list.files("./../../../data/input/GSE193895_RAW/", full.names = T)

Each file contains three columns: GeneID – mouse Ensembl gene Id, Length – transcript length and column, named after sample Id and filled with counts for each gene

1

read_tsv(count_files[1]) %>% head %>% kable

## Rows: 53815 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): GeneID
## dbl (2): Length, norm_1_A_F_run2
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

So, let’s read all the files and create full table of counts:

1
2
3
4
5
6

RNA_counts <- lapply(count_files, function(x) read_tsv(x)) %>% 
              purrr::reduce(dplyr::full_join, by = 'GeneID') %>% 
              rowwise() %>%
              mutate(Glength = max(c_across(starts_with("Length")))) %>% 
              dplyr::select(!starts_with("Length")) %>% 
              ungroup()

View the first few rows of the table:

1

RNA_counts %>% head %>% kable

Fortunately, this data table does not contain any missing values. All samples were processed at the same time and we do not have to deal with them.

1

length(RNA_counts[is.na(RNA_counts)])

## [1] 0

We need this table for all kinds of subsequent analyzes. So, we should write it on the disk:

1
2
3
4

write.table(RNA_counts, 
            file = "./../../../data/input/data//RNA_counts.tsv", 
            row.names = F, 
            sep = "\t")

Add annotation to our counts DT

The GeneID column of our data contains version of each gene. We need gene index without version. So, let’s create new index column:

1
2
3
4

RNA_counts <- RNA_counts %>% 
  mutate(gene_id = sapply(stringr::str_split(RNA_counts$GeneID, 
                                             pattern = "\\.", n = 2),  
                          function(.)  .[1]))

Now we need annotations. AnnotationHub provides interface for all most useful data annotation systems.

1

hub <- AnnotationHub()

## snapshotDate(): 2022-10-31

We could retrieve information about species and data providers:

1

unique(hub$species) %>% head %>% kable

1

hub$species[grep("Mus", hub$species)] %>% head

## [1] "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus" "Mus musculus"
## [6] "Mus musculus"

1
2
3

u <- unique(hub$dataprovider)
o <- order(unique(hub$dataprovider))
u[o] %>% head %>% kable

Obtain list of tables for Mus musculus from Ensembl:

1

query(hub, c("Mus musculus", "Ensembl"))

## AnnotationHub with 2485 records
## # snapshotDate(): 2022-10-31
## # $dataprovider: Ensembl, BioMart, UCSC, FANTOM5,DLRP,IUPHAR,HPRD,STRING,SWI...
## # $species: Mus musculus, mus musculus wsbeij, mus musculus pwkphj, mus musc...
## # $rdataclass: TwoBitFile, GRanges, EnsDb, data.frame, SQLiteFile, list, OrgDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
## #   rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH6066"]]' 
## 
##              title                                              
##   AH6066   | Ensembl Genes                                      
##   AH6093   | Ensembl Genes                                      
##   AH6165   | Ensembl Genes                                      
##   AH6204   | Ensembl Genes                                      
##   AH7567   | Mus_musculus.GRCm38.70.gtf                         
##   ...        ...                                                
##   AH106502 | Mus_musculus_wsbeij.WSB_EiJ_v1.dna_sm.toplevel.2bit
##   AH106503 | Mus_musculus_wsbeij.WSB_EiJ_v1.ncrna.2bit          
##   AH107060 | org.Mm.eg.db.sqlite                                
##   AH107223 | LRBaseDb for Mus musculus (Mouse, v004)            
##   AH109367 | Ensembl 108 EnsDb for Mus musculus

We’ll use the latest one – AH109367 | Ensembl 108 EnsDb for Mus musculus:

1

ensdb <- hub[["AH109367"]]

## loading from cache

Obtain metadata from EnsemblDB for each gene:

1
2

gns <- genes(ensdb)
gns %>% head %>% kable

Add to our data table information about gene and transcript length from Encembl.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# get transcripts
# transcriptsBy includes the UTR and exon regions
txs <- transcriptsBy(ensdb)

# max transcript widths (RNA + UTRs)
# width() is a getter for IRanges objects (<IRanges>)
txwid <- sapply(width(txs), max)

# full DNA gene widths, with names
gnwid <- setNames(width(gns), names(gns))

# txwid and gnwid have different orders
# create table of "TxLen", "GeneLen", "ens_id"
ens.annot <- data.frame(TxLen = txwid, GeneLen = gnwid[names(txwid)]) %>% 
  rownames_to_column(var = "gene_id") 

ens.annot %>% head %>% kable

Merge tables:

1

RNA_counts_ann <- left_join(RNA_counts, ens.annot,  by = "gene_id")

Because TxLen contains not only cDNA sequence it is not useful for us now.

1
2
3

RNA_counts_ann %>% dplyr::select(gene_id, Glength, TxLen, GeneLen) %>% 
  head %>% 
  kable

Add rest of metadata to our table:

1
2

RNA_counts_ann <- left_join(RNA_counts_ann, as.data.frame(gns), by = "gene_id") %>% 
                  unnest(entrezid)

Let’s view annotations of our table:

1
2
3
4

RNA_counts_ann %>% 
  dplyr::select(1,37:53) %>% 
  head %>% 
  kable

So, this is the time to write our table on disk for subsequent uses.

1
2
3
4

write.table(RNA_counts_ann, 
            file = "./../../../data/input/data//RNA_counts_ann.tsv", 
            row.names = F, 
            sep = "\t")

Afterword

In this article we’ve found out how to obtain and prepare RNA seq data. In the next post we will create table of conditions and perform analysis of differential gene expression.

References: