TCGA: gene IDs not appearing and other concerns

Question

I want to download RNA-seq datasets of 4 or 5 different types of cancer from the TCGA to investigate what is happening with my gene of interest. The problem is that I'm a physician and I don't know people used to deal with Big Data to help me.

I'm using a Bioconductor program called TCGAbiolinks and I'm able to download some datasets using the following script:

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("TCGAbiolinks")

browseVignettes("TCGAbiolinks")

library(TCGAbiolinks)
library(dplyr)
library(DT)

query <- GDCquery(project = "TCGA-BRCA", 
                  data.category = "Gene expression",
                  file.type  = "normalized_results",
                  data.type = "Gene expression quantification",
                  experimental.strategy = "RNA-Seq", 
                  barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
                  legacy = TRUE)

GDCdownload(query, method = "api", files.per.chunk = 10)

# GDCprepare transforms data in R format

# The default GDCprepare output is a SummarizedExperiment (SE) object,
# but you can also set to output the data frame.

data <- GDCprepare(query, summarizedExperiment = FALSE)

The problems is that: 1) the gene list is not appearing in the first column; 2) I don't know how many datasets I'll be able to download and analyse without a server.. I'm working with my Macbook..

Can someone familiar with the TCGA website or with the TCGAbiolinks help me with this project? I've been trying to progress in the last months, but the TCGA is a complex website and I'm feeling lost. I use R, I don't know SQL or other programming languages..

Can you show us what head(data) gives you? There ought to be some identifier for the gene - probably an ENST transcript identifier. Datasets for 4-5 types of cancer won't be a problem BTW - they'll take up very little of your storage and if you use data.table (and load as few datasets as you can afford to into memory), not much memory either. How much RAM does your Macbook have? — Ram RS, Commented Aug 17, 2020 at 15:10
thank you very much, @RamRS. Please see the table above, now included in my question. My Mac has 8 GB 2400 MHz DDR4. — FPS, Commented Aug 17, 2020 at 16:01
Try requesting without the summarizedExperiment=FALSE argument - that should give you more metadata. — Ram RS, Commented Aug 17, 2020 at 18:08
it gives me many S4 objects, but I don't know how to work with them.. I don't even know how to open a S4 object.. — FPS, Commented Aug 18, 2020 at 12:15
Take a look at the SummarizedExperiment documentation. You can use assays(data) to get the numbers matrix, and rowData(data) and colData(data) to get information on the features and samples respectively. — Ram RS, Commented Aug 18, 2020 at 19:29

Ram RS · Accepted Answer · 2020-08-19 17:42:00Z

I had to play around a little to understand the object, but here's what I got:

R version 4.0.2 (2020-06-22) -- "Taking Off Again"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TCGAbiolinks)
> query <- GDCquery(project = "TCGA-BRCA",
+                   data.category = "Gene expression",
+                   file.type  = "normalized_results",
+                   data.type = "Gene expression quantification",
+                   experimental.strategy = "RNA-Seq",
+                   barcode = c("TCGA-AN-A046", "TCGA-AC-A23H"),
+                   legacy = TRUE)
--------------------------------------
o GDCquery: Searching in GDC database
--------------------------------------
Genome of reference: hg19
--------------------------------------------
oo Accessing GDC. This might take a while...
--------------------------------------------
ooo Project: TCGA-BRCA
--------------------
oo Filtering results
--------------------
ooo By experimental.strategy
ooo By data.type
ooo By file.type
ooo By barcode
----------------
oo Checking data
----------------
ooo Check if there are duplicated cases
ooo Check if there results for the query
-------------------
o Preparing output
-------------------

> GDCdownload(query, method = "api", files.per.chunk = 10)
Downloading data for project TCGA-BRCA
GDCdownload will download 3 files. A total of 1.306357 MB
Downloading chunk 1 of 1 (3 files, size = 1.306357 MB) as Wed_Aug_19_12_20_11_2020_0.tar.gz

Downloading: 580 kB

> data <- GDCprepare(query, summarizedExperiment = TRUE)
------------------
oo Reading 3 files
------------------
|====================================================|100%                      Completed after 0 s
------------------
oo Merging 3 files
------------------
Accessing grch37.ensembl.org to get gene information
Downloading genome information (try:0) Using: Human genes (GRCh37.p13)
Starting to add information to samples
 => Add clinical information to samples
 => Adding TCGA molecular information from marker papers
 => Information will have prefix 'paper_'
brca subtype information from:doi.org/10.1016/j.ccell.2018.03.014

> head(data)
class: RangedSummarizedExperiment
dim: 6 3
metadata(1): data_release
assays(1): normalized_count
rownames(6): A1BG A2M ... RP11-986E7.7 AADAC
rowData names(3): gene_id entrezgene ensembl_gene_id
colnames(3): TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07
  TCGA-AC-A23H-01A-11R-A157-07
colData names(83): barcode patient ... paper_PARADIGM Clusters
  paper_Pan-Gyn Clusters

I see from this representation that I can use metadata(data), assays(data), etc to get more information.

> head(SummarizedExperiment::assays(data))
List of length 1
names(1): normalized_count

Alright, so this is a named list. I know how to deal with lists: I can use either the index ([[1]]) or the name ($normalized_count)

> head(SummarizedExperiment::assays(data)$normalized_count)
             TCGA-AN-A046-01A-21R-A034-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-01A-11R-A157-07
A1BG                             121.9879                      38.5795                      13.6874
A2M                             9703.8622                   37640.8373                    4767.0490
NAT1                            8488.0727                     114.8582                     270.4753
NAT2                               6.0583                       4.2938                       2.3804
RP11-986E7.7                   33016.6641                   11308.3460                    3894.9639
AADAC                              1.1359                     219.6976                       1.1902

I can now get more information on the samples using colData(). I'm truncating the output below because there are 83 columns.

> SummarizedExperiment::colData(data)
DataFrame with 3 rows and 83 columns
                                                  barcode      patient           sample shortLetterCode          definition sample_submitter_id         sample_type       state
                                              <character>  <character>      <character>     <character>         <character>         <character>         <character> <character>
TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046-01A-21R-A034-07 TCGA-AN-A046 TCGA-AN-A046-01A              TP Primary solid Tumor    TCGA-AN-A046-01A       Primary Tumor    released
TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H-11A-12R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-11A              NT Solid Tissue Normal    TCGA-AC-A23H-11A Solid Tissue Normal    released
TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H-01A-11R-A157-07 TCGA-AC-A23H TCGA-AC-A23H-01A              TP Primary solid Tumor    TCGA-AC-A23H-01A       Primary Tumor    released
                                                        sample_id   is_ffpe sample_type_id  tissue_type submitter_id days_to_collection oct_embedded                pathology_report_uuid
                                                      <character> <logical>    <character>  <character>  <character>          <integer>  <character>                          <character>
TCGA-AN-A046-01A-21R-A034-07 57aff09f-0e97-4e65-a227-3dc7a8516367     FALSE             01 Not Reported TCGA-AN-A046                 27         true 1304FB17-A20A-4EBC-9CC8-1554808AC1F6
TCGA-AC-A23H-11A-12R-A157-07 7df59ca8-7e51-4581-9d7f-8bba0395ce17     FALSE             11 Not Reported TCGA-AC-A23H                478        false                                   NA
TCGA-AC-A23H-01A-11R-A157-07 d7e3b628-d5fd-4e79-9c4a-6409330fb8a7     FALSE             01 Not Reported TCGA-AC-A23H                478        false A7C7D409-D086-4A9B-8C8F-E7E231D5891D

I can get more information on the genes using rowData():

> SummarizedExperiment::rowData(data)
DataFrame with 19947 rows and 3 columns
                  gene_id entrezgene ensembl_gene_id
              <character>  <integer>     <character>
A1BG                 A1BG          1 ENSG00000121410
A2M                   A2M          2 ENSG00000175899
NAT1                 NAT1          9 ENSG00000171428
NAT2                 NAT2         10 ENSG00000156006
RP11-986E7.7 RP11-986E7.7         12 ENSG00000273259
...                   ...        ...             ...
RASAL2-AS1     RASAL2-AS1  100302401 ENSG00000224687
LINC00882       LINC00882  100302640 ENSG00000242759
FTX                   FTX  100302692 ENSG00000230590
TICAM2             TICAM2  100302736 ENSG00000243414
SLC25A5-AS1   SLC25A5-AS1  100303728 ENSG00000224281

You want the assay data as a data.matrix, I guess. That'd be:

> data_as_matrix <- data.matrix(SummarizedExperiment::assays(data)$normalized_count)

My session Info:

> sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] TCGAbiolinks_2.16.1

loaded via a namespace (and not attached):
  [1] bitops_1.0-6                matrixStats_0.56.0          bit64_0.9-7.1               doParallel_1.0.15           RColorBrewer_1.1-2          progress_1.2.2
  [7] httr_1.4.2                  GenomeInfoDb_1.24.2         backports_1.1.8             tools_4.0.2                 R6_2.4.1                    DBI_1.1.0
 [13] BiocGenerics_0.34.0         colorspace_1.4-1            tidyselect_1.1.0            gridExtra_2.3               prettyunits_1.1.1           bit_1.1-15.2
 [19] curl_4.3                    compiler_4.0.2              rvest_0.3.6                 Biobase_2.48.0              xml2_1.3.2                  DelayedArray_0.14.0
 [25] rtracklayer_1.48.0          scales_1.1.1                survMisc_0.5.5              readr_1.3.1                 genefilter_1.70.0           askpass_1.1
 [31] rappdirs_0.3.1              stringr_1.4.0               digest_0.6.25               Rsamtools_2.4.0             foreign_0.8-80              R.utils_2.9.2
 [37] rio_0.5.16                  XVector_0.28.0              pkgconfig_2.0.3             dbplyr_1.4.4                readxl_1.3.1                rlang_0.4.7
 [43] ggthemes_4.2.0              RSQLite_2.2.0               generics_0.0.2              zoo_1.8-8                   jsonlite_1.7.0              BiocParallel_1.22.0
 [49] zip_2.0.4                   car_3.0-8                   dplyr_1.0.0                 R.oo_1.23.0                 RCurl_1.98-1.2              magrittr_1.5
 [55] GenomeInfoDbData_1.2.3      Matrix_1.2-18               Rcpp_1.0.5                  munsell_0.5.0               S4Vectors_0.26.1            abind_1.4-5
 [61] lifecycle_0.2.0             R.methodsS3_1.8.0           stringi_1.4.6               carData_3.0-4               SummarizedExperiment_1.18.1 zlibbioc_1.34.0
 [67] plyr_1.8.6                  BiocFileCache_1.12.0        grid_4.0.2                  blob_1.2.1                  parallel_4.0.2              ggrepel_0.8.2
 [73] forcats_0.5.0               crayon_1.3.4                survminer_0.4.8             lattice_0.20-41             haven_2.3.1                 Biostrings_2.56.0
 [79] splines_4.0.2               GenomicFeatures_1.40.0      annotate_1.66.0             hms_0.5.3                   knitr_1.29                  pillar_1.4.6
 [85] ggpubr_0.4.0                GenomicRanges_1.40.0        ggsignif_0.6.0              codetools_0.2-16            biomaRt_2.44.1              stats4_4.0.2
 [91] XML_3.99-0.5                glue_1.4.1                  downloader_0.4              data.table_1.13.0           vctrs_0.3.2                 selectr_0.4-2
 [97] foreach_1.5.0               cellranger_1.1.0            gtable_0.3.0                openssl_1.4.2               purrr_0.3.4                 tidyr_1.1.0
[103] km.ci_0.5-2                 assertthat_0.2.1            ggplot2_3.3.2               openxlsx_4.1.5              xfun_0.16                   xtable_1.8-4
[109] broom_0.7.0                 rstatix_0.6.0               survival_3.1-12             tibble_3.0.3                iterators_1.0.12            KMsurv_0.1-5
[115] GenomicAlignments_1.24.0    AnnotationDbi_1.50.1        memoise_1.1.0               IRanges_2.22.2              ellipsis_0.3.1

Stack Exchange Network

TCGA: gene IDs not appearing and other concerns

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
r
rna-seq
or ask your own question.

Hot Network Questions

TCGA: gene IDs not appearing and other concerns

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged rrna-seq or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
r
rna-seq
or ask your own question.