RNA-seq workflow: gene-level differential expression rna-seq pipeline 从原始数据到差异分析一条龙

RNA-seq workflow: gene-level exploratory analysis and differential expression

Abstract

Here we walk through an end-to-end gene-level RNA-seq differential expression workflow using Bioconductor packages. We will start from the FASTQ files, show how these were quantified to the reference transcripts, and prepare gene-level count datasets for downstream analysis. We will perform exploratory data analysis (EDA) for quality assessment and to explore the relationship between samples, perform differential gene expression analysis, and visually explore the results.

Contents

• 1 Introduction
  • 1.1 Experimental data
• 2 Preparing quantification input to DESeq2
  • 2.1 Transcript quantification and tximport / tximeta
  • 2.2 Quantifying with Salmon
  • 2.3 Reading in data with tximeta
  • 2.4 DESeq2 import functions
  • 2.5 SummarizedExperiment
  • 2.6 Branching point
• 3 The DESeqDataSet object, sample information and the design formula
  • 3.1 Starting from SummarizedExperiment
  • 3.2 Starting from count matrices
• 4 Exploratory analysis and visualization
  • 4.1 Pre-filtering the dataset
  • 4.2 The variance stabilizing transformation and the rlog
  • 4.3 Sample distances
  • 4.4 PCA plot
  • 4.5 PCA plot using Generalized PCA
  • 4.6 MDS plot
• 5 Differential expression analysis
  • 5.1 Running the differential expression pipeline
  • 5.2 Building the results table
  • 5.3 Other comparisons
  • 5.4 Multiple testing
• 6 Plotting results
  • 6.1 Counts plot
  • 6.2 MA-plot
  • 6.3 Gene clustering
  • 6.4 Independent filtering
  • 6.5 Independent Hypothesis Weighting
• 7 Annotating and exporting results
  • 7.1 Exporting results
  • 7.2 Plotting fold changes in genomic space
• 8 Removing hidden batch effects
  • 8.1 Using SVA with DESeq2
  • 8.2 Using RUV with DESeq2
• 9 Time course experiments
• 10 Session information
• References

R version: R version 4.2.0 RC (2022-04-19 r82224)

Bioconductor version: 3.15

Package: 1.20.0

1.表达矩阵如何获取

可以从geo下载 最后制作成为SummarizedExperiment格式，作为deseq2的输入input

DESeq2 import functions

While the above section described use of Salmon and tximeta, there are many possible inputs to DESeq2, each of which have their own dedicated import functions. The following tools can be used generate or compile count data for use with DESeq2: tximport (Soneson, Love, and Robinson 2015), tximeta (Love et al. 2020), htseq-count (Anders, Pyl, and Huber 2015), featureCounts (Liao, Smyth, and Shi 2014), summarizeOverlaps (Lawrence et al. 2013).

We will next describe the class of object created by tximeta which was saved above as se and gse, and how to create a DESeqDataSet object from this for use with DESeq2 (the other functions above also create a DESeqDataSet).

2.5SummarizedExperiment

The component parts of a SummarizedExperiment object. The assay (pink block) contains the matrix of counts, the rowRanges (blue block) contains information about the genomic ranges and the colData (green block) contains information about the samples. The highlighted line in each block represents the first row (note that the first row of colData lines up with the first column of the assay).

The SummarizedExperiment container is diagrammed in the Figure above and discussed in the latest Bioconductor paper (Huber et al. 2015). In our case, tximeta has created an object gse with three matrices: “counts” - the estimated fragment counts for each gene and sample, “abundance” - the estimated transcript abundances in TPM, and “length” - the effective gene lengths which include changes in length due to biases as well as due to transcript usage. The names of the assays can be examined with assayNames, and the assays themselves are stored as assays (a list of matrices). The first matrix in the list can be pulled out via assay. The rowRanges for our object is the GRanges of the genes (from the left-most position of all the transcripts to the right-most position of all the transcripts). The component parts of the SummarizedExperiment are accessed with an R function of the same name: assay (or assays), rowRanges and colData.

We now will load the full count matrix corresponding to all samples and all data, which is provided in the airway package, and will continue the analysis with the full data object. We can investigate this SummarizedExperiment object by looking at the matrices in the assays slot, the phenotypic data about the samples in colData slot, and the data about the genes in the rowRanges slot.

data(gse)
gse

## class: RangedSummarizedExperiment 
## dim: 58294 8 
## metadata(6): tximetaInfo quantInfo ... txomeInfo txdbInfo
## assays(3): counts abundance length
## rownames(58294): ENSG00000000003.14 ENSG00000000005.5 ...
##   ENSG00000285993.1 ENSG00000285994.1
## rowData names(1): gene_id
## colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
## colData names(3): names donor condition

The counts are the first matrix, so we can examine them with just assay:

assayNames(gse)

## [1] "counts"    "abundance" "length"

head(assay(gse), 3)

##                    SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14    708.164    467.962    900.992    424.368   1188.295
## ENSG00000000005.5       0.000      0.000      0.000      0.000      0.000
## ENSG00000000419.12    455.000    510.000    604.000    352.000    583.000
##                    SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14   1090.668    805.929    599.337
## ENSG00000000005.5       0.000      0.000      0.000
## ENSG00000000419.12    773.999    409.999    499.000

colSums(assay(gse))

## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 
##   21100805   19298584   26145537   15688246   25268618   31891456   19683767 
## SRR1039521 
##   21813903

The rowRanges, when printed, shows the ranges for the first five and last five genes:

rowRanges(gse)

## GRanges object with 58294 ranges and 1 metadata column:
##                      seqnames              ranges strand |            gene_id
##                         <Rle>           <IRanges>  <Rle> |        <character>
##   ENSG00000000003.14     chrX 100627109-100639991      - | ENSG00000000003.14
##    ENSG00000000005.5     chrX 100584802-100599885      + |  ENSG00000000005.5
##   ENSG00000000419.12    chr20   50934867-50958555      - | ENSG00000000419.12
##   ENSG00000000457.13     chr1 169849631-169894267      - | ENSG00000000457.13
##   ENSG00000000460.16     chr1 169662007-169854080      + | ENSG00000000460.16
##                  ...      ...                 ...    ... .                ...
##    ENSG00000285990.1    chr14   19244904-19269380      - |  ENSG00000285990.1
##    ENSG00000285991.1     chr6 149817937-149896011      - |  ENSG00000285991.1
##    ENSG00000285992.1     chr8   47129262-47132628      + |  ENSG00000285992.1
##    ENSG00000285993.1    chr18   46409197-46410645      - |  ENSG00000285993.1
##    ENSG00000285994.1    chr10   12563151-12567351      + |  ENSG00000285994.1
##   -------
##   seqinfo: 25 sequences (1 circular) from hg38 genome

The rowRanges also contains metadata about the sequences (chromosomes in our case) in the seqinfo slot:

seqinfo(rowRanges(gse))

## Seqinfo object with 25 sequences (1 circular) from hg38 genome:
##   seqnames seqlengths isCircular genome
##   chr1      248956422      FALSE   hg38
##   chr2      242193529      FALSE   hg38
##   chr3      198295559      FALSE   hg38
##   chr4      190214555      FALSE   hg38
##   chr5      181538259      FALSE   hg38
##   ...             ...        ...    ...
##   chr21      46709983      FALSE   hg38
##   chr22      50818468      FALSE   hg38
##   chrX      156040895      FALSE   hg38
##   chrY       57227415      FALSE   hg38
##   chrM          16569       TRUE   hg38

The colData for the SummarizedExperiment reflects the data.frame that was provided to the tximeta function for importing the quantification data. Here we can see that there are columns indicating sample names, as well as the donor ID, and the treatment condition (treated with dexamethasone or untreated).

colData(gse)

## DataFrame with 8 rows and 3 columns
##                 names    donor     condition
##              <factor> <factor>      <factor>
## SRR1039508 SRR1039508  N61311  Untreated    
## SRR1039509 SRR1039509  N61311  Dexamethasone
## SRR1039512 SRR1039512  N052611 Untreated    
## SRR1039513 SRR1039513  N052611 Dexamethasone
## SRR1039516 SRR1039516  N080611 Untreated    
## SRR1039517 SRR1039517  N080611 Dexamethasone
## SRR1039520 SRR1039520  N061011 Untreated    
## SRR1039521 SRR1039521  N061011 Dexamethasone

2.6Branching point

At this point, we have counted the fragments which overlap the genes in the gene model we specified. This is a branching point where we could use a variety of Bioconductor packages for exploration and differential expression of the count data, including edgeR (Robinson, McCarthy, and Smyth 2009), limma with the voom method (Law et al. 2014), DSS (Wu, Wang, and Wu 2013), EBSeq (Leng et al. 2013) and baySeq (Hardcastle and Kelly 2010). Schurch et al. (2016) compared performance of different statistical methods for RNA-seq using a large number of biological replicates and can help users to decide which tools make sense to use, and how many biological replicates are necessary to obtain a certain sensitivity. We will continue using DESeq2 (Love, Huber, and Anders 2014). The SummarizedExperiment object is all we need to start our analysis. In the following section we will show how to use it to create the data object used by DESeq2.

In DESeq2, the custom class is called DESeqDataSet. It is built on top of the SummarizedExperiment class, and it is easy to convert SummarizedExperiment objects into DESeqDataSet objects, which we show below. One of the two main differences is that the assay slot is instead accessed using the counts accessor function, and the DESeqDataSet class enforces that the values in this matrix are non-negative integers.

A second difference is that the DESeqDataSet has an associated design formula. The experimental design is specified at the beginning of the analysis, as it will inform many of the DESeq2 functions how to treat the samples in the analysis (one exception is the size factor estimation, i.e., the adjustment for differing library sizes, which does not depend on the design formula). The design formula tells which columns in the sample information table (colData) specify the experimental design and how these factors should be used in the analysis.

First, let’s examine the columns of the colData of gse. We can see each of the columns just using the $ directly on the SummarizedExperiment or DESeqDataSet.

gse$donor

## [1] N61311  N61311  N052611 N052611 N080611 N080611 N061011 N061011
## Levels: N052611 N061011 N080611 N61311

gse$condition

## [1] Untreated     Dexamethasone Untreated     Dexamethasone Untreated    
## [6] Dexamethasone Untreated     Dexamethasone
## Levels: Untreated Dexamethasone

We can rename our variables if we want. Let’s use cell to denote the donor cell line, and dex to denote the treatment condition.

gse$cell <- gse$donor
gse$dex <- gse$condition

We can also change the names of the levels. It is critical when one renames levels to not change the order. Here we will rename "Untreated" as "untrt" and "Dexamethasone" as "trt":

levels(gse$dex)

## [1] "Untreated"     "Dexamethasone"

# when renaming levels, the order must be preserved!
levels(gse$dex) <- c("untrt", "trt")

The simplest design formula for differential expression would be ~ condition, where condition is a column in colData(dds) that specifies which of two (or more groups) the samples belong to. For the airway experiment, we will specify ~ cell + dex meaning that we want to test for the effect of dexamethasone (dex) controlling for the effect of different cell line (cell).

Note: it is prefered in R that the first level of a factor be the reference level (e.g. control, or untreated samples). In this case, when the colData table was assembled the untreated samples were already set as the reference, but if this were not the case we could use relevel as shown below. While levels(...) <- above was simply for renaming the character strings associated with levels, relevel is a very different function, which decides how the variables will be coded, and how contrasts will be computed. For a two-group comparison, the use of relevel to change the reference level would flip the sign of a coefficient associated with a contrast between the two groups.

library("magrittr")
gse$dex %<>% relevel("untrt")
gse$dex

## [1] untrt trt   untrt trt   untrt trt   untrt trt  
## Levels: untrt trt

%<>% is the compound assignment pipe-operator from the magrittr package, the above line of code is a concise way of saying:

gse$dex <- relevel(gse$dex, "untrt")

For running DESeq2 models, you can use R’s formula notation to express any fixed-effects experimental design. Note that DESeq2 uses the same formula notation as, for instance, the lm function of base R. If the research aim is to determine for which genes the effect of treatment is different across groups, then interaction terms can be included and tested using a design such as ~ group + treatment + group:treatment. See the manual page for ?results for more examples. We will show how to use an interaction term to test for condition-specific changes over time in a time course example below.

In the following sections, we will demonstrate the construction of the DESeqDataSet from two starting points:  构造输入数据input 

GSE82158_Bulk-seq_bleomycin_IM_AM_macrophage_ontogeny\\GSE82158_macrophage_ontogeny_norm_counts‘_YoungLeeyou的博客-CSDN博客

• from a SummarizedExperiment object
• from a count matrix and a sample information table

For a full example of using the HTSeq Python package for read counting, please see the pasilla vignette. For an example of generating the DESeqDataSet from files produced by htseq-count, please see the DESeq2 vignette.

Starting from SummarizedExperiment

Again, we can quickly check the millions of fragments that could be mapped by Salmon to the genes (the second argument of round tells how many decimal points to keep).

round( colSums(assay(gse)) / 1e6, 1 )

## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 
##       21.1       19.3       26.1       15.7       25.3       31.9       19.7 
## SRR1039521 
##       21.8

Once we have our fully annotated SummarizedExperiment object, we can construct a DESeqDataSet object from it that will then form the starting point of the analysis. We add an appropriate design for the analysis:

library("DESeq2")

dds <- DESeqDataSet(gse, design = ~ cell + dex)

3.2Starting from count matrices

In this section, we will show how to build an DESeqDataSet supposing we only have a count matrix and a table of sample information.

Note: if you have prepared a SummarizedExperiment you should skip this section. While the previous section would be used to construct a DESeqDataSet from a SummarizedExperiment, here we first extract the individual object (count matrix and sample info) from the SummarizedExperiment in order to build it back up into a new object – only for demonstration purposes. In practice, the count matrix would either be read in from a file or perhaps generated by an R function like featureCounts from the Rsubread package (Liao, Smyth, and Shi 2014).

The information in a SummarizedExperiment object can be accessed with accessor functions. For example, to see the actual data, i.e., here, the fragment counts, we use the assay function. (The head function restricts the output to the first few lines.)

countdata <- round(assays(gse)[["counts"]])
head(countdata, 3)

##                    SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003.14        708        468        901        424       1188
## ENSG00000000005.5           0          0          0          0          0
## ENSG00000000419.12        455        510        604        352        583
##                    SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003.14       1091        806        599
## ENSG00000000005.5           0          0          0
## ENSG00000000419.12        774        410        499

In this count matrix, each row represents a gene, each column a sequenced RNA library, and the values give the estimated counts of fragments that were probabilistically assigned to the respective gene in each library by Salmon. We also have information on each of the samples (the columns of the count matrix). If you’ve imported the count data in some other way, for example loading a pre-computed count matrix, it is very important to check manually that the columns of the count matrix correspond to the rows of the sample information table.

coldata <- colData(gse)

We now have all the ingredients to prepare our data object in a form that is suitable for analysis, namely:

• countdata: a table with the fragment counts
• coldata: a table with information about the samples

To now construct the DESeqDataSet object from the matrix of counts and the sample information table, we use:

ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = coldata,
                                 design = ~ cell + dex)

We will continue with the object generated from the SummarizedExperiment section.