3.1.1 Sample types

A key feature distinguishing the pancancer-atlas project from TCGA is the availability of data from normal tissue or metastatic or recurrent tumor samples. Codes are used to distinguish the different sources:

BiocOncoTK::pancan_sampTypeMap

##   SampleTypeLetterCode                                      SampleType
## 1                  TAM                           Additional Metastatic
## 2                  TAP                        Additional - New Primary
## 3                   TR                           Recurrent Solid Tumor
## 4                   TB Primary Blood Derived Cancer - Peripheral Blood
## 5                   TM                                      Metastatic
## 6                   NT                             Solid Tissue Normal
## 7                   TP                             Primary solid Tumor

3.1.2 SummarizedExperiments per assay per tumor or other sample type

The following code will run if you have a valid setting for environment variable CGC_BILLING, to allow BiocOncoTK::pancan_BQ() to generate a proper BigQueryConnection.

library(BiocOncoTK)
if (nchar(Sys.getenv("CGC_BILLING"))>0) {
 pcbq = pancan_BQ() # basic connection
 BRCA_mir = restfulSE::pancan_SE(pcbq)
}

The result is

> BRCA_mir
class: SummarizedExperiment 
dim: 743 1068 
metadata(0):
assays(1): assay
rownames(743): hsa-miR-30d-3p hsa-miR-486-3p ... hsa-miR-525-3p
  hsa-miR-892b
rowData names(0):
colnames(1068): TCGA-LD-A7W6 TCGA-BH-A18I ... TCGA-E9-A1N9 TCGA-B6-A0X0
colData names(746): bcr_patient_uuid bcr_patient_barcode ...
  bilirubin_upper_limit days_to_last_known_alive

3.1.3 Subsetting to normal

To shift attention to the normal tissue samples provided, use

 BRCA_mir_nor = restfulSE::pancan_SE(pcbq, assaySampleTypeCode="NT")

to find

class: SummarizedExperiment 
dim: 743 90 
metadata(0):
assays(1): assay
rownames(743): hsa-miR-7641 hsa-miR-135a-5p ... hsa-miR-1323
  hsa-miR-520d-5p
rowData names(0):
colnames(90): TCGA-BH-A18P TCGA-BH-A18S ... TCGA-E9-A1N6 TCGA-E9-A1N9
colData names(746): bcr_patient_uuid bcr_patient_barcode ...
  bilirubin_upper_limit days_to_last_known_alive

The intersection of the colnames from the two SummarizedExperiments thus formed (patients contributing both solid tumor and matched normal) has length 89.

3.1.4 Shifting focus to another tissue/assay

You need to know what type of sample has been assayed for the tumor type of interest.

Here is how you find the candidates.

bqcon %>% tbl(pancan_longname("rnaseq")) %>% filter(Study=="GBM") %>% 
   group_by(SampleTypeLetterCode) %>% summarise(n=n())

To get RNA-seq on recurrent GBM samples:

pancan_SE(bqcon, colDFilterValue="GBM", tumorFieldValue="GBM", 
  assayDataTableName=pancan_longname("rnaseq"), 
  assaySampleTypeCode="TR", assayFeatureName="Symbol", 
  assayValueFieldName="normalized_count")

3.1.5 Multiassay experiments per tumor

Suppose we want to work with the mRNA, RPPA, 27k/450k merged methylation and miRNA data together. We can invoke pancan_SE again, specifying the appropriate tables and fields.

BRCA_mrna = pancan_SE(pcbq,
   assayDataTableName = pancan_longname("rnaseq"),
   assayFeatureName = "Entrez",
   assayValueFieldName = "normalized_count")
BRCA_rppa = pancan_SE(pcbq,
   assayDataTableName = pancan_longname("RPPA"),
   assayFeatureName = "Protein",
   assayValueFieldName = "Value")
BRCA_meth = pancan_SE(pcbq,
   assayDataTableName = pancan_longname("27k")[2],
   assayFeatureName = "ID",
   assayValueFieldName = "Beta")

After obtaining the clinical data for BRCA with

library(dplyr)
library(magrittr)
clinBRCA = pcbq %>% tbl(pancan_longname("clinical")) %>% 
  filter(acronym=="BRCA") %>% as.data.frame() 
rownames(clinBRCA) = clinBRCA[,2]
clinDF = DataFrame(clinBRCA)

we use

library(MultiAssayExperiment)
brcaMAE = MultiAssayExperiment(
  ExperimentList(rnaseq=BRCA_mrna, meth=BRCA_meth, rppa=BRCA_rppa,
    mirna=BRCA_mir),colData=clinDF)

to generate brcaMAE. No assay data are present in this object, but data are retrieved on request.

> brcaMAE
A MultiAssayExperiment object of 4 listed
 experiments with user-defined names and respective classes. 
 Containing an ExperimentList class object of length 4: 
 [1] rnaseq: SummarizedExperiment with 20531 rows and 1097 columns 
 [2] meth: SummarizedExperiment with 22601 rows and 1067 columns 
 [3] rppa: SummarizedExperiment with 259 rows and 873 columns 
 [4] mirna: SummarizedExperiment with 743 rows and 1068 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

It is convenient to check for sample availability for the different assays using upsetSamples in MultiAssayExperiment.

The upset diagram for brcaMAE, showing sample availability per assay.

3.1.6 Reassessing immune infiltrate signatures in relation to microsatellite instability

The following code produces figure 1 of the restfulSE supplement.

library(BiocOncoTK)

 infilGenes = c(`PD-L1`="CD274", `PD-L2`="PDCD1LG2", CD8A="CD8A")
 tumcodes = c("COAD", "STAD", "UCEC")
 combs = expand.grid(tumcode=tumcodes, ali=names(infilGenes),
    stringsAsFactors=FALSE)
 combs$sym = infilGenes[combs$ali]

 bq = pancan_BQ()
 exprByMSI = function(bq, tumcode, genesym, alias) {
  print(tumcode)
  if (missing(alias)) alias=genesym
  ex = bindMSI(buildPancanSE(bq, tumcode, assay="RNASeqv2"))
  ex = replaceRownames(ex)
  data.frame(
   patient_barcode=colnames(ex),
   acronym=tumcode,
   symbol = genesym,
   alias = alias,
   log2ex=log2(as.numeric(SummarizedExperiment::assay(ex[genesym,]))+1),
   msicode = ifelse(ex$msiTest >= 4, ">=4", "<4"))
 }
 allshow = lapply(1:nrow(combs), function(x) exprByMSI(bq, combs$tumcode[x],
    combs$sym[x], combs$ali[x]))
 
 rr = do.call(rbind, allshow)
 
 library(ggplot2)
 png(file="microsatpan2.png")
 ggplot(rr,
    aes(msicode, log2ex)) + geom_boxplot() +
    facet_grid(acronym~alias) +
    ylab("log2(normalized expr. + 1)") +
    xlab("microsatellite instability score")
 dev.off()

3.5.1 Mutation data

First we get an overview of the content:

muttab = con %>% tbl("Mutation_calls")
length(muttab %>% colnames())
muttab %>% select(Cell_line_primary_name, Hugo_Symbol, 
   Variant_Classification, cDNA_Change)%>% glimpse()

## [1] 53

Now let’s filter by NRAS and get a feel for how many observations are returned per cell line.

nrastab = muttab %>% select(Variant_Classification, Hugo_Symbol, 
    Cell_line_primary_name, CCLE_name) %>%
     filter(Hugo_Symbol == "NRAS") %>% group_by(Hugo_Symbol) 
nrastab %>% summarise(n=n())
nrasdf = nrastab %>% as.data.frame()

We need to carve up the CCLE name to get the organ.

spl = function(x) {
  z = strsplit(x, "_")
  fir = vapply(z, function(x)x[1], character(1))
  rest = vapply(z, function(x) paste(x[-1], collapse="_"), character(1))
  list(fir, rest)
}
nrasdf$organ = spl(nrasdf$CCLE_name)[[2]]

head(nrasdf)

##   Variant_Classification Hugo_Symbol Cell_line_primary_name
## 1      Missense_Mutation        NRAS                SNU-387
## 2      Missense_Mutation        NRAS                SNU-719
## 3      Missense_Mutation        NRAS                 SNU-81
## 4      Missense_Mutation        NRAS                SW 1271
## 5      Missense_Mutation        NRAS                   TF-1
## 6      Missense_Mutation        NRAS                  THP-1
##                                 CCLE_name                              organ
## 1                            SNU387_LIVER                              LIVER
## 2                          SNU719_STOMACH                            STOMACH
## 3                   SNU81_LARGE_INTESTINE                    LARGE_INTESTINE
## 4                             SW1271_LUNG                               LUNG
## 5  TF1_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE HAEMATOPOIETIC_AND_LYMPHOID_TISSUE
## 6 THP1_HAEMATOPOIETIC_AND_LYMPHOID_TISSUE HAEMATOPOIETIC_AND_LYMPHOID_TISSUE

table(nrasdf$organ)

## 
##                  AUTONOMIC_GANGLIA             CENTRAL_NERVOUS_SYSTEM 
##                                  2                                  1 
##                        ENDOMETRIUM HAEMATOPOIETIC_AND_LYMPHOID_TISSUE 
##                                  3                                 25 
##                    LARGE_INTESTINE                              LIVER 
##                                  2                                  2 
##                               LUNG                              OVARY 
##                                  9                                  2 
##                               SKIN                        SOFT_TISSUE 
##                                  6                                  2 
##                            STOMACH                            THYROID 
##                                  1                                  1 
##          UPPER_AERODIGESTIVE_TRACT                      URINARY_TRACT 
##                                  2                                  3

prim_names = as.character(nrasdf$Cell_line_primary_name)

3.5.2 Expression data

Let’s obtain the expression of AHR for these NRAS-mutated cell lines.

ccexp = con %>% tbl("AffyU133_RMA_expression") 
ccexp %>% glimpse()
ccexp %>% select(Cell_line_primary_name, RMA_normalized_expression,
    HGNC_gene_symbol) %>% filter(HGNC_gene_symbol == "AHR") %>% 
    filter(Cell_line_primary_name %in% nrasdf$Cell_line_primary_name) %>%
    as.data.frame() -> NRAS_AHR
head(NRAS_AHR)

##   Cell_line_primary_name RMA_normalized_expression HGNC_gene_symbol
## 1                SNU-719                 11.859660              AHR
## 2                 SNU-81                 10.684210              AHR
## 3              SK-MEL-30                  9.901173              AHR
## 4               KU-19-19                  8.772602              AHR
## 5                  L-363                  3.704958              AHR
## 6                 DND-41                  4.110147              AHR

3.5.3 Drug responsiveness data from CCLE, using pogos

The pogos package (submitted, see github.com/vjcitn/pogos) includes software to query pharmacodb.pmgenomics.ca. We will use this to develop drug-response profiles for PD-0325901.

library(pogos)
ccleNRAS = DRTraceSet(NRAS_AHR[,1], drug="PD-0325901")
plot(ccleNRAS)

We’ll define a responsiveness method, that takes a function f that is applied to the responses component of the dose-response profile.

responsiveness = function (x, f) 
{
    r = sapply(slot(x, "traces"), function(x) f(slot(slot(x,"DRProfiles")[[1]],"responses")))
    data.frame(Cell_line_primary_name = slot(x,"cell_lines"), resp = r, 
        drug = slot(x,"drug"), dataset = x@dataset)
}

The activity area for a compound in this design is defined as

AA = function(x) sum((pmax(0, x/100)))
head(rr <- responsiveness(ccleNRAS, AA))

##   Cell_line_primary_name   resp       drug dataset
## 1              SK-MEL-30 2.8200 PD-0325901    CCLE
## 2                  L-363 8.0061 PD-0325901    CCLE
## 3                SK-N-AS 4.0700 PD-0325901    CCLE
## 4                SNU-387 7.5157 PD-0325901    CCLE
## 5                   Mino 6.8010 PD-0325901    CCLE
## 6                 MOLP-8 5.4010 PD-0325901    CCLE

summary(rr$resp)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.920   3.730   5.178   4.958   5.927   8.006

This is based on the supplement to Barretina et al. 2012. (There a slightly different formula in the addendum which uses notation that includes multiplying by a factor of i for dose index level i.)

Let’s merge the responsiveness data with the expression data for gene AHR.

rexp = merge(rr, NRAS_AHR)
rexp[1:2,]

##   Cell_line_primary_name   resp       drug dataset RMA_normalized_expression
## 1                  639-V 5.5290 PD-0325901    CCLE                  4.711943
## 2                    697 6.1672 PD-0325901    CCLE                  4.803876
##   HGNC_gene_symbol
## 1              AHR
## 2              AHR

3.5.4 CLUE

The CLUE platform is an interface to results of work on the connectivity map at Broad Institute. Usage of functions in this toolkit requires an API key, which can be acquired through registration at clue.io. Set the environment variable CLUE_KEY so that it can be found by Sys.getenv to use default key parameter to functions described here.

A basic purpose of the interface to CLUE is to allow identification of gene signatures of perturbations in specific cellular contexts.

We have serialized data on cell lines and perturbagens available in the GSE70138 snapshot of LINCS.

data(cell_70138)
names(cell_70138)

##  [1] "cell_id"                 "cell_type"              
##  [3] "base_cell_id"            "precursor_cell_id"      
##  [5] "modification"            "sample_type"            
##  [7] "primary_site"            "subtype"                
##  [9] "original_growth_pattern" "provider_catalog_id"    
## [11] "original_source_vendor"  "donor_age"              
## [13] "donor_sex"               "donor_ethnicity"

table(cell_70138$primary_site)

## 
##                               -666                            adipose 
##                                  4                                  2 
##                  autonomic ganglia                              blood 
##                                  1                                  1 
##                               bone                             breast 
##                                  4                                  9 
##             central nervous system                        endometrium 
##                                  3                                  3 
## haematopoietic and lymphoid tissue                             kidney 
##                                  8                                  6 
##                    large intestine                              liver 
##                                 17                                  3 
##                               lung                             muscle 
##                                 13                                  3 
##                              ovary                           pancreas 
##                                  6                                  2 
##                           prostate                               skin 
##                                  5                                  6 
##                            stomach                    vascular system 
##                                  1                                  1

data(pert_70138)
dim(pert_70138)

## [1] 2170    3

names(pert_70138)

## [1] "pert_id"    "pert_iname" "pert_type"

A number of API services have demonstration query expressions available in the package:

cd = clueDemos()
names(cd)

## [1] "cells"         "genes"         "profiles"      "perts"        
## [5] "sigs"          "rep_drug_moas" "pcls"

cd$sigs

## $where
## $where$pert_desc
## [1] "imatinib"
## 
## $where$cell_id
## [1] "MCF7"

We use query_clue to query a service. Here we ask for perturbagens that have EGFR among their targets. We’ll retrieve a single ‘gold’ signature identifier.

if (nchar(Sys.getenv("CLUE_KEY"))>0) {
lkbytarg = query_clue(service="perts", filter=list(where=list(target="EGFR")))
print(names(lkbytarg[[1]]))
sig1 = lkbytarg[[1]]$sig_id_gold[1]
}

Now we obtain the metadata about this signature.

if (nchar(Sys.getenv("CLUE_KEY"))>0) {
sig1d = query_clue(service="sigs", filter=list(where=list(sig_id=sig1)))
print(names(sig1d[[1]]))
print(head(sig1d[[1]]$pert_iname)) # perturbagen
print(head(sig1d[[1]]$cell_id))  # cell type
print(head(sig1d[[1]]$dn50_lm))  # some downregulated genes among the landmark
print(head(sig1d[[1]]$up50_lm))  # some upregulated genes among the landmark
}

# use pertClasses() to get names of perturbagen classes in Clue
if (nchar(Sys.getenv("CLUE_KEY"))>0) {
tuinh = query_clue("perts", 
   filter=list(where=list(pcl_membership=list(inq=list("CP_HDAC_INHIBITOR"))))) 
inames_tu = sapply(tuinh, function(x)x$pert_iname)

npcSigs = query_clue(service="sigs", filter=list(where=list(cell_id="NPC")))
length(npcSigs)
gns = lapply(npcSigs, function(x) x$up50_lm)
perts = lapply(npcSigs, function(x) x$pert_iname)
touse = which(perts %in% inames_tu)
rec = names(tab <- sort(table(unlist(gns[touse])),decreasing=TRUE)[1:5])
cbind(select(org.Hs.eg.db, keys=rec, columns="SYMBOL"), n=as.numeric(tab))
}