sparseAssay

The most straightforward matrix representation of a RaggedExperiment will return a matrix with the number of rows equal to the total number of ranges defined across all samples. i.e. the 5 rows of the sparseAssay result:

sparseAssay(ragexp)

##   sample1 sample2
## A       3      NA
## B       4      NA
## C       5      NA
## D      NA      11
## E      NA      12

correspond to the ranges of the unlisted GRangesList:

unlist(grl)

## GRanges object with 5 ranges and 2 metadata columns:
##             seqnames    ranges strand |     score        type
##                <Rle> <IRanges>  <Rle> | <integer> <character>
##   sample1.A     chr1      1-10      - |         3    germline
##   sample1.B     chr1      8-14      + |         4     somatic
##   sample1.C     chr1     15-18      + |         5    germline
##   sample2.D     chr1      1-10      - |        11    germline
##   sample2.E     chr1     11-18      + |        12     somatic
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

The rownames of the sparseAssay result are equal to the names of the GRanges elements. The values in the matrix returned by sparseAssay correspond to the first columns of the mcols of each GRangesList element, in this case the “score” column.

Note, this is the default assay() method of RaggedExperiment:

assay(ragexp, "score")

##   sample1 sample2
## A       3      NA
## B       4      NA
## C       5      NA
## D      NA      11
## E      NA      12

assay(ragexp, "type")

##   sample1    sample2   
## A "germline" NA        
## B "somatic"  NA        
## C "germline" NA        
## D NA         "germline"
## E NA         "somatic"

compactAssay

The dimensions of the compactAssay result differ from that of the sparseAssay result only if there are identical ranges in different samples. Identical ranges are placed in the same row in the output. Ranges with any difference in start, end, or strand, will be kept on different rows. Non-disjoint ranges are not collapsed.

compactAssay(ragexp)

##              sample1 sample2
## chr1:8-14:+        4      NA
## chr1:11-18:+      NA      12
## chr1:15-18:+       5      NA
## chr1:1-10:-        3      11

compactAssay(ragexp, "type")

##              sample1    sample2   
## chr1:8-14:+  "somatic"  NA        
## chr1:11-18:+ NA         "somatic" 
## chr1:15-18:+ "germline" NA        
## chr1:1-10:-  "germline" "germline"

Note that row names are constructed from the ranges, and the names of the GRanges vectors are lost, unlike in the sparseAssay result.

disjoinAssay

This function is similar to compactAssay except the rows are disjoint¹ ranges. Elements of the matrix are summarized by applying the simplifyDisjoin functional argument to assay values of overlapping ranges.

disjoinAssay(ragexp, simplifyDisjoin = mean)

##              sample1 sample2
## chr1:8-10:+        4      NA
## chr1:11-14:+       4      12
## chr1:15-18:+       5      12
## chr1:1-10:-        3      11

qreduceAssay

The qreduceAssay function is the most complicated but likely the most useful of the RaggedExperiment Assay functions. It requires you to provide a query argument that is a GRanges vector, and the rows of the resulting matrix correspond to the elements of this GRanges. The returned matrix will have dimensions length(query) by ncol(x). Elements of the resulting matrix correspond to the overlap of the i th query range in the j th sample, summarized according to the simplifyReduce functional argument. This can be useful, for example, to calculate per-gene copy number or mutation status by providing the genomic ranges of every gene as the query.

The simplifyReduce argument in qreduceAssay allows the user to summarize overlapping regions with a custom method for the given “query” region of interest. We provide one for calculating a weighted average score per query range, where the weight is proportional to the overlap widths between overlapping ranges and a query range.

Note that there are three arguments to this function. See the documentation for additional details.

weightedmean <- function(scores, ranges, qranges)
{
    isects <- pintersect(ranges, qranges)
    sum(scores * width(isects)) / sum(width(isects))
}

The call to qreduceAssay calculates the overlaps between the ranges of each sample:

grl

## GRangesList object of length 2:
## $sample1
## GRanges object with 3 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        type
##        <Rle> <IRanges>  <Rle> | <integer> <character>
##   A     chr1      1-10      - |         3    germline
##   B     chr1      8-14      + |         4     somatic
##   C     chr1     15-18      + |         5    germline
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
## 
## $sample2
## GRanges object with 2 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        type
##        <Rle> <IRanges>  <Rle> | <integer> <character>
##   D     chr1      1-10      - |        11    germline
##   E     chr1     11-18      + |        12     somatic
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

with the query ranges (an arbitrary set is defined here for demonstration): First create a demonstration “query” region of interest:

(query <- GRanges(c("chr1:1-14:-", "chr1:15-18:+")))

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1      1-14      -
##   [2]     chr1     15-18      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

using the simplifyReduce function to resolve overlapping ranges and return a matrix with rows corresponding to the query:

qreduceAssay(ragexp, query, simplifyReduce = weightedmean)

##              sample1 sample2
## chr1:1-14:-        3      11
## chr1:15-18:+       5      12

Input data format

For demonstration, we consider CNV calls as obtained with PennCNV from SNP-chip data in a Brazilian cattle breed (da Silva et al., 2016).

Here, we use a data subset and only consider CNV calls on chromosome 1 and 2, for which there are roughly 3000 CNV calls as obtained for 711 samples.

data.dir <- system.file("extdata", package="CNVRanger")
call.file <- file.path(data.dir, "Silva16_PONE_CNV_calls.csv")
calls <- read.csv(call.file, as.is=TRUE)
nrow(calls)

## [1] 3000

head(calls)

##    chr start    end    NE_id state
## 1 chr1 16947  45013 NE001423     3
## 2 chr1 36337  67130 NE001426     3
## 3 chr1 16947  36337 NE001428     3
## 4 chr1 36337 105963 NE001519     3
## 5 chr1 36337  83412 NE001534     3
## 6 chr1 36337  83412 NE001648     3

length(unique(calls[,"NE_id"]))

## [1] 711

Representation as a GRangesList

We group the calls by sample ID, resulting in a GRangesList. Each element of the list corresponds to a sample, and contains the genomic coordinates of the CNV calls for this sample (along with the copy number state in the State metadata column).

grl <- makeGRangesListFromDataFrame(calls, 
    split.field="NE_id", keep.extra.columns=TRUE)
grl

## GRangesList object of length 711:
## $NE001357
## GRanges object with 5 ranges and 1 metadata column:
##       seqnames            ranges strand |     state
##          <Rle>         <IRanges>  <Rle> | <integer>
##   [1]     chr1   4569526-4577608      * |         3
##   [2]     chr1 15984544-15996851      * |         1
##   [3]     chr1 38306432-38330161      * |         3
##   [4]     chr1 93730576-93819471      * |         0
##   [5]     chr2 40529044-40540747      * |         3
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $NE001358
## GRanges object with 1 range and 1 metadata column:
##       seqnames              ranges strand |     state
##          <Rle>           <IRanges>  <Rle> | <integer>
##   [1]     chr1 105042452-105233446      * |         1
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $NE001359
## GRanges object with 4 ranges and 1 metadata column:
##       seqnames            ranges strand |     state
##          <Rle>         <IRanges>  <Rle> | <integer>
##   [1]     chr1   4569526-4577608      * |         3
##   [2]     chr1 31686841-31695808      * |         0
##   [3]     chr1 93730576-93819471      * |         0
##   [4]     chr2   2527718-2535261      * |         0
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## ...
## <708 more elements>

The advantage of representing the CNV calls as a GRangesList is that it allows to leverage the comprehensive set of operations on genomic regions implemented in the GenomicRanges package - for instance, sorting of the calls according to their genomic coordinates.

grl <- sort(grl)
grl

## GRangesList object of length 711:
## $NE001357
## GRanges object with 5 ranges and 1 metadata column:
##       seqnames            ranges strand |     state
##          <Rle>         <IRanges>  <Rle> | <integer>
##   [1]     chr1   4569526-4577608      * |         3
##   [2]     chr1 15984544-15996851      * |         1
##   [3]     chr1 38306432-38330161      * |         3
##   [4]     chr1 93730576-93819471      * |         0
##   [5]     chr2 40529044-40540747      * |         3
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $NE001358
## GRanges object with 1 range and 1 metadata column:
##       seqnames              ranges strand |     state
##          <Rle>           <IRanges>  <Rle> | <integer>
##   [1]     chr1 105042452-105233446      * |         1
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## $NE001359
## GRanges object with 4 ranges and 1 metadata column:
##       seqnames            ranges strand |     state
##          <Rle>         <IRanges>  <Rle> | <integer>
##   [1]     chr1   4569526-4577608      * |         3
##   [2]     chr1 31686841-31695808      * |         0
##   [3]     chr1 93730576-93819471      * |         0
##   [4]     chr2   2527718-2535261      * |         0
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths
## 
## ...
## <708 more elements>

Representation as a RaggedExperiment

An alternative matrix-like representation of the CNV calls can be obtained with the RaggedExperiment data class. It resembles in many aspects the SummarizedExperiment data class for storing gene expression data as e.g. obtained with RNA-seq.

ra <- RaggedExperiment(grl)
ra

## class: RaggedExperiment 
## dim: 3000 711 
## assays(1): state
## rownames: NULL
## colnames(711): NE001357 NE001358 ... NE003967 NE003968
## colData names(0):

As apparent from the dim slot of the object, it stores the CNV calls in the rows and the samples in the columns. Note that the CN state is now represented as an assay matrix which can be easily accessed and subsetted.

assay(ra[1:5,1:5])

##                        NE001357 NE001358 NE001359 NE001360 NE001361
## chr1:4569526-4577608          3       NA       NA       NA       NA
## chr1:15984544-15996851        1       NA       NA       NA       NA
## chr1:38306432-38330161        3       NA       NA       NA       NA
## chr1:93730576-93819471        0       NA       NA       NA       NA
## chr2:40529044-40540747        3       NA       NA       NA       NA

As with SummarizedExperiment objects, additional information for the samples are annotated in the colData slot. For example, we annotate the steer weight and its feed conversion ratio (FCR) using simulated data. Feed conversion ratio is the ratio of dry matter intake to live-weight gain. A typical range of feed conversion ratios is 4.5-7.5 with a lower number being more desirable as it would indicate that a steer required less feed per pound of gain.

weight <- rnorm(ncol(ra), mean=1100, sd=100)
fcr <- rnorm(ncol(ra), mean=6, sd=1.5)
colData(ra)$weight <- round(weight, digits=2)
colData(ra)$fcr <- round(fcr, digits=2)
colData(ra)

## DataFrame with 711 rows and 2 columns
##             weight       fcr
##          <numeric> <numeric>
## NE001357   1146.49      5.98
## NE001358   1178.38      7.76
## NE001359   1072.19      5.74
## NE001360   1043.34      7.33
## NE001361   1166.39      5.21
## ...            ...       ...
## NE003962    983.84      5.87
## NE003963     983.1      6.06
## NE003966   1103.54      6.75
## NE003967   1334.75      6.25
## NE003968   1099.83      5.87

Trimming low-density areas

Here, we use the approach from CNVRuler to summarize CNV calls to CNV regions (see Figure 1 in Kim et al., 2012 for an illustration of the approach). This trims low-density areas as defined by the density argument, which is set here to <10% of the number of calls within a summarized region.

cnvrs <- populationRanges(grl, density=0.1)
cnvrs

## GRanges object with 303 ranges and 2 metadata columns:
##         seqnames              ranges strand |      freq        type
##            <Rle>           <IRanges>  <Rle> | <numeric> <character>
##     [1]     chr1        16947-111645      * |       103        gain
##     [2]     chr1     1419261-1630187      * |        18        gain
##     [3]     chr1     1896112-2004603      * |       218        gain
##     [4]     chr1     4139727-4203274      * |         1        gain
##     [5]     chr1     4554832-4577608      * |        23        gain
##     ...      ...                 ...    ... .       ...         ...
##   [299]     chr2 136310067-136322489      * |         2        loss
##   [300]     chr2 136375337-136386940      * |         1        loss
##   [301]     chr2 136455546-136466040      * |         1        loss
##   [302]     chr2 136749793-136802493      * |        39        both
##   [303]     chr2 139194749-139665914      * |        58        both
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Note that CNV frequency (number of samples overlapping each region) and CNV type (gain, loss, or both) have also been annotated in the columns freq and type, respectively.

Identifying recurrent regions

In particular in cancer, it is important to distinguish driver from passenger mutations, i.e. to distinguish meaningful events from random background aberrations. The GISTIC method identifies those regions of the genome that are aberrant more often than would be expected by chance, with greater weight given to high amplitude events (high-level copy-number gains or homozygous deletions) that are less likely to represent random aberrations (Beroukhim et al., 2007).

By setting est.recur=TRUE, we deploy a GISTIC-like significance estimation

cnvrs <- populationRanges(grl, density=0.1, est.recur=TRUE)
cnvrs

## GRanges object with 303 ranges and 3 metadata columns:
##         seqnames              ranges strand |      freq        type
##            <Rle>           <IRanges>  <Rle> | <numeric> <character>
##     [1]     chr1        16947-111645      * |       103        gain
##     [2]     chr1     1419261-1630187      * |        18        gain
##     [3]     chr1     1896112-2004603      * |       218        gain
##     [4]     chr1     4139727-4203274      * |         1        gain
##     [5]     chr1     4554832-4577608      * |        23        gain
##     ...      ...                 ...    ... .       ...         ...
##   [299]     chr2 136310067-136322489      * |         2        loss
##   [300]     chr2 136375337-136386940      * |         1        loss
##   [301]     chr2 136455546-136466040      * |         1        loss
##   [302]     chr2 136749793-136802493      * |        39        both
##   [303]     chr2 139194749-139665914      * |        58        both
##                      pvalue
##                   <numeric>
##     [1] 0.00980392156862742
##     [2]   0.107843137254902
##     [3]                   0
##     [4]   0.558823529411765
##     [5]  0.0882352941176471
##     ...                 ...
##   [299]   0.236111111111111
##   [300]   0.421296296296296
##   [301]   0.421296296296296
##   [302]  0.0588235294117647
##   [303]  0.0392156862745098
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

and filter for recurrent CNVs that exceed a significance threshold of 0.05.

cnvrs[cnvrs$pvalue < 0.05]

## GRanges object with 17 ranges and 3 metadata columns:
##        seqnames              ranges strand |      freq        type
##           <Rle>           <IRanges>  <Rle> | <numeric> <character>
##    [1]     chr1        16947-111645      * |       103        gain
##    [2]     chr1     1896112-2004603      * |       218        gain
##    [3]     chr1   15984544-15996851      * |       116        loss
##    [4]     chr1   31686841-31695808      * |       274        loss
##    [5]     chr1   69205418-69219486      * |        46        loss
##    ...      ...                 ...    ... .       ...         ...
##   [13]     chr2   97882695-97896935      * |        80        loss
##   [14]     chr2 124330343-124398570      * |        39        loss
##   [15]     chr2 135096060-135271140      * |        84        gain
##   [16]     chr2 135290754-135553033      * |        83        gain
##   [17]     chr2 139194749-139665914      * |        58        both
##                     pvalue
##                  <numeric>
##    [1] 0.00980392156862742
##    [2]                   0
##    [3]  0.0185185185185185
##    [4] 0.00462962962962965
##    [5]  0.0416666666666666
##    ...                 ...
##   [13]  0.0231481481481481
##   [14]  0.0462962962962963
##   [15]  0.0196078431372549
##   [16]  0.0294117647058824
##   [17]  0.0392156862745098
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Working with individual CNV and RNA-seq assays

We demonstrate the functionality by loading RNA-seq read count data from skeletal muscle samples for 183 Nelore cattle steers, which we analyzed for CNV-expression effects as previously described (Geistlinger et al., 2018).

rseq.file <- file.path(data.dir, "counts_cleaned.txt")
rcounts <- read.delim(rseq.file, row.names=1, stringsAsFactors=FALSE)
rcounts <- as.matrix(rcounts)
dim(rcounts)

## [1] 939 183

rcounts[1:5, 1:5]

##                    NE001407 NE001408 NE001424 NE001439 NE001448
## ENSBTAG00000000088       64       65      233      135      134
## ENSBTAG00000000160       20       28       50       13       18
## ENSBTAG00000000176      279      373      679      223      417
## ENSBTAG00000000201      252      271      544      155      334
## ENSBTAG00000000210      268      379      486      172      443

For demonstration, we restrict analysis to the 939 genes on chromosome 1 and 2, and store the RNA-seq expression data in a SummarizedExperiment.

library(SummarizedExperiment)
rse <- SummarizedExperiment(assays=list(rcounts=rcounts), 
                rowRanges=granges(sel.genes)[rownames(rcounts)])
rse

## class: RangedSummarizedExperiment 
## dim: 939 183 
## metadata(0):
## assays(1): rcounts
## rownames(939): ENSBTAG00000000088 ENSBTAG00000000160 ...
##   ENSBTAG00000048210 ENSBTAG00000048228
## rowData names(0):
## colnames(183): NE001407 NE001408 ... NE003840 NE003843
## colData names(0):

Assuming distinct modes of action, effects observed in the CNV-expression analysis are typically divided into (i) local effects (cis), where expression changes coincide with CNVs in the respective genes, and (ii) distal effects (trans), where CNVs supposedly affect trans-acting regulators such as transcription factors.

However, due to power considerations and to avoid detection of spurious effects, stringent filtering of (i) not sufficiently expressed genes, and (ii) CNV regions with insufficient sample size in groups deviating from 2n, should be carried out when testing for distal effects. Local effects have a clear spatial indication and the number of genes locating in or close to a CNV region of interest is typically small; testing for differential expression between CN states is thus generally better powered for local effects and less stringent filter criteria can be applied.

In the following, we carry out CNV-expression association analysis by providing the CNV regions to test (cnvrs), the individual CNV calls (grl) to determine per-sample CN state in each CNV region, the RNA-seq read counts (rse), and the size of the genomic window around each CNV region (window). The window argument thereby determines which genes are considered for testing for each CNV region and is set here to 1 Mbp.

Further, use the min.cpm and min.samples arguments to exclude from the analysis (i) genes with fewer than min.cpm reads per million reads mapped (cpm, counts per million), and (ii) CNV regions with fewer than min.samples samples in a group deviating from 2n.

res <- cnvEQTL(cnvrs, grl, rse, window="1Mbp", verbose=TRUE)

## Restricting analysis to 179 intersecting samples

## Preprocessing RNA-seq data ...

## Summarizing per-sample CN state in each CNV region

## Excluding 286 cnvrs not satisfying min.samples threshold

## Analyzing 12 regions with >=1 gene in the given window

## 1 of 12

## 2 of 12

## 3 of 12

## 4 of 12

## 5 of 12

## 6 of 12

## 7 of 12

## 8 of 12

## 9 of 12

## 10 of 12

## 11 of 12

## 12 of 12

res

## DataFrame with 89 rows and 5 columns
##                         CNVR               Gene                 logFC
##                  <character>        <character>             <numeric>
## 1          chr1:16947-111645 ENSBTAG00000018278     0.194859708204693
## 2          chr1:16947-111645 ENSBTAG00000021997   -0.0812496004002579
## 3          chr1:16947-111645 ENSBTAG00000020035    0.0745028563584966
## 4          chr1:16947-111645 ENSBTAG00000011528    0.0118392582964246
## 5          chr1:16947-111645 ENSBTAG00000012594  -0.00661188778069866
## ...                      ...                ...                   ...
## 85  chr2:135290754-135553033 ENSBTAG00000030340   -0.0310183382350007
## 86  chr2:135290754-135553033 ENSBTAG00000008314   -0.0148596021181995
## 87  chr2:135290754-135553033 ENSBTAG00000008309    0.0172997593964661
## 88  chr2:135290754-135553033 ENSBTAG00000003403  -0.00336281264078722
## 89  chr2:135290754-135553033 ENSBTAG00000013282 -0.000223673376098313
##                  PValue         AdjPValue
##               <numeric>         <numeric>
## 1   0.00731514635939691 0.400665393143022
## 2     0.179616021014321 0.841564657598637
## 3     0.740325113667423 0.997718543512711
## 4     0.915565410373996 0.997718543512711
## 5     0.943057146391629 0.997718543512711
## ...                 ...               ...
## 85    0.850536851523179 0.997718543512711
## 86    0.869368789923888 0.997718543512711
## 87    0.885514362762372 0.997718543512711
## 88    0.987950415154231 0.997718543512711
## 89    0.997718543512711 0.997718543512711

The resulting DataFrame contains in the first column the CNV regions tested. The second column contains the genes tested in the genomic window around each CNV region, and the other columns report (i) log2 fold change with respect to the 2n group, (ii) edgeR’s DE p-value, and (iii) the (per default) Benjamini-Hochberg adjusted p-value.

Working with TCGA data stored in a MultiAssayExperiment

In the previous section, we individually prepared the CNV and RNA-seq data for CNV-expression association analysis. In the following, we demonstrate how to perform an integrated preparation of the two assays when stored in a MultiAssayExperiment We therefore consider glioblastoma GBM data from The Cancer Genome Atlas TCGA, which can conveniently be accessed with the curatedTCGAData package.

library(curatedTCGAData)
suppressMessages(
    gbm <- curatedTCGAData("GBM", 
        assays=c("GISTIC_Peaks", "CNVSNP", "RNASeq2GeneNorm"), 
        dry.run=FALSE)
)
gbm

## A MultiAssayExperiment object of 3 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 3: 
##  [1] GBM_CNVSNP-20160128: RaggedExperiment with 146852 rows and 1104 columns 
##  [2] GBM_RNASeq2GeneNorm-20160128: SummarizedExperiment with 20501 rows and 166 columns 
##  [3] GBM_GISTIC_Peaks-20160128: RangedSummarizedExperiment with 68 rows and 577 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

The returned MultiAssayExperiment contains three assays:

• the SNP-based CNV calls stored in a RaggedExperiment (GBM_CNVSNP),
• the recurrent CNV regions summarized across the population using the GISTIC method (GBM_GISTIC_Peaks), and
• the normalized RNA-seq gene expression values in a SummarizedExperiment (GBM_RNASeq2GeneNorm).

To annotate the genomic coordinates of the genes measured in the RNA-seq assay, we use the function symbolsToRanges from the TCGAutils package. For demonstration, we restrict the analysis to chromosome 1 and 2.

library(TCGAutils)
gbm <- suppressMessages(symbolsToRanges(gbm, unmapped=FALSE))
for(i in 1:3) 
{
    genome(rowRanges(gbm[[i]])) <- "hg19"
    seqlevelsStyle(rowRanges(gbm[[i]])) <- "UCSC"
    ind <- as.character(seqnames(rowRanges(gbm[[i]]))) %in% c("chr1", "chr2")
    gbm[[i]] <- gbm[[i]][ind,]
}
gbm

## A MultiAssayExperiment object of 3 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 3: 
##  [1] GBM_CNVSNP-20160128: RaggedExperiment with 17818 rows and 1104 columns 
##  [2] GBM_GISTIC_Peaks-20160128: RangedSummarizedExperiment with 12 rows and 577 columns 
##  [3] GBM_RNASeq2GeneNorm-20160128_ranged: RangedSummarizedExperiment with 2939 rows and 166 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

We now restrict the analysis to intersecting patients of the three assays using MultiAssayExperiment’s intersectColumns function, and select Primary Solid Tumor samples using the splitAssays function from the TCGAutils package.

gbm <- intersectColumns(gbm)
sampleTables(gbm)

## $`GBM_CNVSNP-20160128`
## 
##  01  02  10  11 
## 154  13 146   1 
## 
## $`GBM_GISTIC_Peaks-20160128`
## 
##  01 
## 154 
## 
## $`GBM_RNASeq2GeneNorm-20160128_ranged`
## 
##  01  02 
## 147  13

data(sampleTypes)
sampleTypes

##    Code                                        Definition
## 1    01                               Primary Solid Tumor
## 2    02                             Recurrent Solid Tumor
## 3    03   Primary Blood Derived Cancer - Peripheral Blood
## 4    04      Recurrent Blood Derived Cancer - Bone Marrow
## 5    05                          Additional - New Primary
## 6    06                                        Metastatic
## 7    07                             Additional Metastatic
## 8    08                        Human Tumor Original Cells
## 9    09        Primary Blood Derived Cancer - Bone Marrow
## 10   10                              Blood Derived Normal
## 11   11                               Solid Tissue Normal
## 12   12                                Buccal Cell Normal
## 13   13                           EBV Immortalized Normal
## 14   14                                Bone Marrow Normal
## 15   15                                    sample type 15
## 16   16                                    sample type 16
## 17   20                                   Control Analyte
## 18   40 Recurrent Blood Derived Cancer - Peripheral Blood
## 19   50                                        Cell Lines
## 20   60                          Primary Xenograft Tissue
## 21   61                Cell Line Derived Xenograft Tissue
## 22   99                                    sample type 99
##    Short.Letter.Code
## 1                 TP
## 2                 TR
## 3                 TB
## 4               TRBM
## 5                TAP
## 6                 TM
## 7                TAM
## 8               THOC
## 9                TBM
## 10                NB
## 11                NT
## 12               NBC
## 13              NEBV
## 14               NBM
## 15              15SH
## 16              16SH
## 17             CELLC
## 18               TRB
## 19              CELL
## 20                XP
## 21               XCL
## 22              99SH

gbm <- splitAssays(gbm, sampleCodes="01")
gbm

## A MultiAssayExperiment object of 3 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 3: 
##  [1] 01_GBM_CNVSNP-20160128: RaggedExperiment with 17818 rows and 154 columns 
##  [2] 01_GBM_GISTIC_Peaks-20160128: RangedSummarizedExperiment with 12 rows and 154 columns 
##  [3] 01_GBM_RNASeq2GeneNorm-20160128_ranged: RangedSummarizedExperiment with 2939 rows and 147 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

We transform CN estimates into CN states where states are encoded as:

• 0: homozygous deletion (2-copy loss)
• 1: heterozygous deletion (1-copy loss)
• 2: normal diploid state
• 3: 1-copy gain
• 4: amplification (>= 2-copy gain)

ind <- grep("CNVSNP", names(gbm))
state <- round(mcols(gbm[[ind]])$Segment_Mean)
state[state > 2] <- 2
state[state < -2] <- -2
mcols(gbm[[ind]])$State <- state + 2
mcols(gbm[[ind]]) <- mcols(gbm[[ind]])[,3:1]
table(mcols(gbm[[ind]])$State)

## 
##     0     1     2     3     4 
##  3270  3025 10223   856   444

The data is now ready for CNV-expression association analysis, where we find only two CNV regions with sufficient sample size for testing using the default value of 10 for the minSamples argument.

res <- cnvEQTL(cnvrs="01_GBM_GISTIC_Peaks-20160128", 
    calls="01_GBM_CNVSNP-20160128", 
    rcounts="01_GBM_RNASeq2GeneNorm-20160128_ranged", 
    data=gbm, window="1Mbp", verbose=TRUE)

## harmonizing input:
##   removing 154 sampleMap rows not in names(experiments)

## Preprocessing RNA-seq data ...

## Summarizing per-sample CN state in each CNV region

## Excluding 10 cnvrs not satisfying min.samples threshold

## Analyzing 2 regions with >=1 gene in the given window

## 1 of 2

## 2 of 2

res

## DataFrame with 39 rows and 5 columns
##                     CNVR        Gene              logFC
##              <character> <character>          <numeric>
## 1   chr1:3394251-6475685       RPL22 -0.725652140752693
## 2   chr1:3394251-6475685      CAMTA1 -0.710092209593344
## 3   chr1:3394251-6475685      KLHL21 -0.804994378329488
## 4   chr1:3394251-6475685     DNAJC11 -0.639331827106962
## 5   chr1:3394251-6475685       PHF13 -0.644451457052169
## ...                  ...         ...                ...
## 35  chr1:7908902-8336254       VAMP3 -0.380592228642377
## 36  chr1:7908902-8336254        H6PD -0.487860462635182
## 37  chr1:7908902-8336254        PER3 -0.496046999162807
## 38  chr1:7908902-8336254      ERRFI1 -0.503783986501104
## 39  chr1:7908902-8336254      SLC2A5 -0.445492206932205
##                   PValue            AdjPValue
##                <numeric>            <numeric>
## 1   9.61335285501219e-09 1.87460380672738e-07
## 2    2.5708616701305e-07  2.1115304553107e-06
## 3   5.49625845092592e-06 3.57256799310185e-05
## 4   6.45920365913183e-06  3.5986991815163e-05
## 5   1.15203955333233e-05 4.99217139777344e-05
## ...                  ...                  ...
## 35   0.00155925324719568  0.00357711039062539
## 36   0.00420206699725692  0.00819403064465099
## 37    0.0280814360824631   0.0456323336340025
## 38    0.0500141322286501   0.0722717568907876
## 39    0.0643968913279903   0.0866027159238491

Prepare input files using PureCN and MuTect

PureCN requires three raw inputs and subsequent preprocessing of them before performing purity and ploidy estimation. Here are the list of three raw input files, followed by the simple workflow diagram.

Five intermediate files

Preprocessing of three input files creates five major intermediate files that are used directly for purity and ploidy estimation with PureCN.R command line tool.

java -jar mutect.jar \
    --analysis_type MuTect \
    -R hg38.fasta \
    --dbsnp $DBSNP_VCF \
    --cosmic $COSMIC_VCF \
    -I:tumor $BAM_TUMOR  \
    -o $OUT/${SAMPLEID}_mutect_stats.txt \
    -vcf $OUT/${SAMPLEID}_mutect.vcf

Rscript $PURECN/PureCN.R \
    --out $PCN_OUT/PureCN/tumor_only/$SAMPLEID \
    --tumor $PCN_OUT/tumor_cov/${SAMPLEID}_coverage_loess.txt \
    --SAMPLEID ${SAMPLEID} \
    --vcf $MUTECT_OUT/stat_tumor_only/${SAMPLEID}_mutect.vcf \
    --statsfile $MUTECT_OUT/stat_tumor_only/${SAMPLEID}_mutect_stats.txt \
    --normaldb $PCN_OUT/normalDB/normalDB_hg38.rds \
    --normal_panel $PCN_OUT/normalDB/mapping_bias_hg38.rds \
    --intervals $INPUT/bed/baits_hg38_intervals.txt \
    --intervalweightfile $PCN_OUT/normalDB/interval_weights_hg38.txt \
    --snpblacklist hg38_simpleRepeats.bed \
    --genome hg38 \
    --force --postoptimize --seed 123

1. Overview

The colors visualize the copy number fitting score from low (blue) to high (red). The numbers indicate the ranks of the local optima.

plotAbs(ret, type="overview")

2. Log-ratio histogram

This figure displays a histogram of tumor vs. normal copy number log-ratios for the maximum likelihood solution (number 1 of ‘Overview’ plot). The height of a bar in this plot is proportional to the fraction of the genome falling into the particular log-ratio copy number range.

plotAbs(ret, 1, type="hist")

## NULL

createCurationFile(file.rds)
read.csv(system.file("extdata/Sample1_PureCN.csv", package = "CNVWorkshop"))

##   Sampleid Purity   Ploidy Sex Contamination Flagged Failed Curated
## 1  Sample1   0.65 1.734572   ?             0    TRUE  FALSE   FALSE
##                          Comment
## 1 EXCESSIVE LOSSES;EXCESSIVE LOH

3. Amplification and deletion

gene.calls = callAlterations(ret)
head(gene.calls)

##            chr     start       end C seg.mean seg.id number.targets
## EIF2A     chr3 150270143 150301699 6   1.3351      5             13
## MIR548H2  chr3 151531951 151538241 6   1.3351      5              3
## GPNMB     chr7  23286477  23313844 7   1.4723     17             11
## SH2D4B   chr10  82298088  82403838 0  -1.5118     22              8
## SLC35G1  chr10  95653791  95661248 0  -1.3002     24              3
## CYP2C9   chr10  96698440  96748786 0  -1.3002     24              9
##           gene.mean   gene.min   gene.max focal breakpoints          type
## EIF2A     1.5159387  0.7544062  2.1392167  TRUE           0 AMPLIFICATION
## MIR548H2  0.8261594  0.4682628  1.2000253  TRUE           0 AMPLIFICATION
## GPNMB     1.4490939  0.9397072  1.8200775  TRUE           0 AMPLIFICATION
## SH2D4B   -1.5141982 -1.9281007 -1.2607230 FALSE           0      DELETION
## SLC35G1  -1.5988526 -1.8066816 -1.4493561 FALSE           0      DELETION
## CYP2C9   -1.1901085 -2.1922992 -0.9875222 FALSE           0      DELETION
##          num.snps  M M.flagged  loh
## EIF2A           0 NA        NA   NA
## MIR548H2        0 NA        NA   NA
## GPNMB           0 NA        NA   NA
## SH2D4B          2  0      TRUE TRUE
## SLC35G1         0 NA        NA   NA
## CYP2C9          0 NA        NA   NA

4. Genomic regions with LOH

loh = callLOH(ret)
head(loh)

##    chr     start       end arm C M                       type   seg.mean
## 1 chr1 114515871 121535434   p 2 0 WHOLE ARM COPY-NEUTRAL LOH  0.1209814
## 2 chr1 124535434 248085104   q 2 0 WHOLE ARM COPY-NEUTRAL LOH  0.1209814
## 3 chr2  10262881  92326171   p 1 0              WHOLE ARM LOH -0.4656416
## 4 chr2  95326171 231775678   q 1 0              WHOLE ARM LOH -0.4656416
## 5 chr2 236403331 239039169   q 2 0           COPY-NEUTRAL LOH  0.0193000
## 6 chr3  11888043  90504854   p 2 1                             0.1052000
##   num.mark num.snps M.flagged maf.expected maf.observed
## 1      902       13     FALSE    0.1750000    0.1756325
## 2      902       13     FALSE    0.1750000    0.1756325
## 3      685       13     FALSE    0.2592593    0.2788990
## 4      685       13     FALSE    0.2592593    0.2788990
## 5       88        2      TRUE    0.1750000    0.1736667
## 6      408        1      TRUE    0.5000000    0.4712710

Example of TCGA-OV and TCGA-LUAD datasets

1. Purity and Ploidy
2. TMB (Tumor mutation burden) : TMB is the number of somatic mutations carried by tumor cells.
3. Mutational signatures : Mutational signatures are characteristic combinations of mutation types arising from specific mutagenesis processes, which can provide information on cause and treatment option for different types of cancer.