Best Way to Download Tcga Data Hoseq
- ane. Introduction
- 2. Licensing and availability
- 3. Installation
- iii.1. Install from CRAN
- 3.2. Install from download package
- three.3. Arrangement compatibility
- 4. Import OMICs profiles
- 4.one. Import factor expression
- 4.1.1. mRNA-microarray gene expression data
- four.1.2. RNA-Seq gene expression data
- 4.2. Import miRNA expression
- 4.2.1. miRNA-array data
- iv.ii.2. miRNA-Seq information
- 4.3. Import mutation
- 4.four. Import methylation expression
- four.5. Import re-create number changes
- 4.one. Import factor expression
- five. Miscellaneous Functions
- 5.one. Combine multiple OMICs profiles
- five.two. Separating OMICs profiles by sample types
- 5.iii. Getting matched tumor-normal OMICs profiles
- 6. Case: Canonical correlation assay
- 7. Appendix
- seven.ane. Appendix A
- 7.2. Appendix B
- seven.3. Appendix C
- eight. Parcel Archive
- 9. Packet In Python
i. Introduction
TCGA2STAT enables users to easily download TCGA data directly into a format ready for statistical analysis in the R surround. The package imports and processes both molecular profiles and clinical information for more than 30 cancer types profiled with dissimilar high-throughput platforms (Appendix A), including microarray, side by side generation sequencing, methylation assortment, SNP array, and assortment-CGH. The package automatically combines molecular profiles and clinical data into a data matrix ready for supervised analysis. Farther, the package provides users a uncomplicated function to merge data from multiple molecular profiles or platforms for integrated statistical assay. About chiefly, all of the parcel functionality is available from ane uncomplicated and user-friendly interface that does not crave whatsoever domain-specific noesis to utilize.
The information imported by this parcel is the version-stamped standardized data sets hosted and maintained by the Broad Plant GDAC Firehose. All imported information will be the latest version available from Firehose; typically, this is the aggregated Level three data from all sample batches and public clinical data equally available from TCGA Data Portal.
This package depends on XML (Lang 2013) to parse the contents and paths while importing data from Firehose. Additionally, it requires CNTools (Zhang) to process the segmented CNV or CNA data into gene-level data.
This parcel relies on the system'south capability to decompress the downloaded data which is a tar. Thus, users running R in Windows system that does not support tar will receive an error message.
2. Licensing and availability
The TCGA2STAT package and the underlying code in the bundle are distributed under the GPL2.0 license. Users of the package are gratis to apply and redistribute the package merely editing is prohibited. The TCGA2STAT packet is bachelor from The Comprehensive R Archive Network (CRAN) and the package's Github.
Usage of this package to import the TCGA data and employ the data in downstream analysis constitutes to agreement to TCGA data usage policy.
3. Installation
three.1. Install from CRAN
The TCGA2STAT package is available from The Comprehensive R Archive Network (CRAN).
iii.2. Install from download bundle
To install TCGA2STAT from the bundle archive file obtained from the packet's Github:
> install.packages("TCGA2STAT_1.0.tar.gz", repos = Null, type = "source") iii.3. System compatibility
This package is intended to be used in Unix and Mac Os, where tar archive utility is supported natively. In Windows, if tar is non supported, users will be given an error message stating Fault: TAR is not installed in the arrangement. Data unzip failed.
To enable tar support in Windows, users could take the following steps:
-
Install Cygwin
-
In R, set the system environment to path for
tarandgzip, for example:
> Sys.setenv(TAR="C:/cygwin64/bin/tar", R_GZIPCMD="C:/cygwin64/bin/gzip") iv. Import OMICs profiles
TCGA2STAT enables users to download any type of data via 1 single, uniform interface: the getTCGA role. The type of data to be downloaded can be specified via three parameters:
-
disease- acronym for the cancer type -
information.type- they blazon of omics-profiles -
type- the specific type of measurement for the omics-profile
For example, the following command volition return the RPKM values of genes profiled using RNA-Sequencing from ovarian cancer patients.
> rnaseq.ov <- getTCGA(disease="OV", information.blazon="RNASeq", type="RPKM") Upon running, a message will be appear indicating that the RNAseq data is being imported. And, upon successful execution of the command, the number of genes imported will be given in another message; for instance:
RNAseq data volition exist imported! This may have some time! 19990 genes have been imported! The returned object rnaseq.ov will be a list containing three elements:
-
dat- The omics-profiles data matrix with genes in rows and patients in columns; in this example, this is a matrix of RPKM values. -
clinical- Clinical data matrix; in this case, Cypher is returned every bit the clinical data is not specified (default). -
merged.dat- If clinical information is imported, a information matrix with both omics-profiles and clinical covariates (overall survival by default); patients are in the rows with the sample ID in the showtime column, followed by clinical covariates, and so the omics-profiles.
Equally an case, the internal structure of the object returned from the to a higher place command is as follows:
> str(rnaseq.ov) List of 3 $ dat : num [i:19990, 1:299] 0.0539 0 0.4204 3.6155 5.3981 ... ..- attr(*, "dimnames")=List of two .. ..$ : chr [one:19990] "AADACL3" "AADACL4" "ABCA4" "ABCB10" ... .. ..$ : chr [1:299] "TCGA-04-1348-01A-01R-1565-thirteen-2" "TCGA-04-1357-01A-01R-1565-thirteen-2" "TCGA-04-1362-01A-01R-1565-thirteen-2" "TCGA-04-1364-01A-01R-1565-13-2" ... $ clinical : Zip $ merged.dat: Nix The information matrix of the omics profiles returned are of dimension factor 10 patients; the row names are gene symbols and the cavalcade names are samples ID (BCR code):
> head(rnaseq.ov$dat[, 1:3]) TCGA-04-1348-01A-01R-1565-13-2 TCGA-04-1357-01A-01R-1565-13-2 TCGA-04-1362-01A-01R-1565-xiii-2 AADACL3 0.0539 0.2341 0.0573 AADACL4 0.0000 0.0000 0.0000 ABCA4 0.4204 0.1647 0.2120 ABCB10 iii.6155 7.1614 3.1549 ABCD3 5.3981 3.8884 seven.2949 ABL2 one.6614 two.4809 one.7997
To obtain default clinical data together with the omics profiles, specify clinical=TRUE. Consider the example beneath:
# Get the RPKM of genes along with all clinical data, and combine the RPKM with overall survival (Bone) > rnaseq_os.ov <- getTCGA(illness="OV", information.type="RNASeq", type="RPKM", clinical=TRUE) # Await at the merged RPKM-Bone data matrix > dim(rnaseq_os.ov$merged.dat) [ane] 295 19993 > head(rnaseq_os.ov$merged.dat[,i:5]) bcr status Bone AADACL3 AADACL4 1 TCGA-04-1348 1 1483 0.0539 0.0000 2 TCGA-04-1357 0 NA 0.2341 0.0000 3 TCGA-04-1362 one 1348 0.0573 0.0000 four TCGA-04-1364 i 1024 0.0000 0.0242 5 TCGA-04-1365 0 2329 0.0177 0.0091 6 TCGA-04-1514 one 1720 0.0019 0.0000 To customize the clinical variables, specify cvars to the desired clinical variable. Consider the case below:
# Become the RPKM gene expression profiles forth with all clinical data, and combine the expression with historic period > rnaseq_age.ov <- getTCGA(disease="OV", data.blazon="RNASeq", type="RPKM", clinical=Truthful, cvars="yearstobirth") # Look at the merged RPKM-age data matrix > head(rnaseq_age.ov$merged.dat[,1:5]) bcr YEARSTOBIRTH AADACL3 AADACL4 ABCA4 ane TCGA-04-1348 44 0.0539 0.0000 0.4204 2 TCGA-04-1357 52 0.2341 0.0000 0.1647 three TCGA-04-1362 59 0.0573 0.0000 0.2120 4 TCGA-04-1364 61 0.0000 0.0242 0.0560 5 TCGA-04-1365 87 0.0177 0.0091 0.0757 six TCGA-04-1514 45 0.0019 0.0000 0.1713 The list of acronyms for all bachelor cancer types is given in Appendix A at the end of this vignette. A summary of all values permitted for both data.type and type is given in Appendix B. A list of values for cvar and their cursory descriptions is given in Appendix C. However, since the clinical variables vary from one cancer type to another, the list in Appendix C is not the complete list. Users interested in a specific clinical covariate for the detail cancer type should check the availability of the clinical covariate from the column-proper name of clinical data matrix in the list returned past getTCGA office.
In the subsections below, we volition give examples of how to download each type of omics information along with a short description of the data obtained.
four.1. Import cistron expression
The ii major platforms used to profile cistron expression in TCGA cancer patients are microarray and RNA-Sequencing. Our package allows users to download data for both platforms if data is bachelor.
iv.1.i. mRNA-microarray factor expression data
Gene expression data profiled by microarray are bachelor in different forms depending on the blazon of microarray employed. Specifically, three forms are available:
- Agilent 244K Custom Gene Expression G4502A-07 (
blazon="G450") - Affymetrix Human Genome U133A 2.0 Array (
type="U133") - Affymetrix Human being Exon ane.0 ST Assortment (
type="Huex")
The information imported from these platforms are level-three gene-level expression information. Please also accept note that but Agilent G450 is available for every cancer type; Affymetrix U133A and Human Exon ST Array are only available for sure cancers.
Agilent 244K Custom Cistron Expression G4502A
# Get Agilent G450 expression for ovarian cancer patients > exp.ov <- getTCGA(disease="OV", data.type="mRNA_Array", type="G450") # Look at the information > dim(exp.ov$dat) [ane] 17814 597 > caput(exp.ov$dat[,1:3]) TCGA-01-0630-11A-01R-0363-07 TCGA-01-0631-11A-01R-0363-07 TCGA-01-0633-11A-01R-0363-07 ELMO2 0.092000 0.0081875 -0.0953125 CREB3L1 0.469500 0.1473333 0.0235000 RPS11 0.965600 0.7915000 1.2695000 PNMA1 0.815400 0.8960000 0.8314000 MMP2 -0.969125 -0.0452500 -0.4230000 C10orf90 -1.887200 -1.9418000 -1.9832000 > This cistron expression information from Agilent G450 should have ~17814 genes. The expression profiles are lowess-normalized log-ratio values that unremarkably range from -10 to 10, as shown in the following histogram (Fig.1):
Fig.1. Histogram of Agilent G450 expression data.
Affymetrix Homo Genome U133A 2.0 Array
# Go Affymetrix Human Genome U133A expression for ovarian cancer patients > u133a.ov <- getTCGA(affliction="OV", data.type="mRNA_Array", type="U133") # Await at the data > dim(u133a.ov$dat) [1] 12042 593 > head(u133a.ov$dat[,1:3]) TCGA-01-0628-11A-01R-0362-01 TCGA-01-0630-11A-01R-0362-01 TCGA-01-0631-11A-01R-0362-01 AACS five.772263 6.428424 five.574237 FSTL1 8.220414 8.496560 8.710241 ELMO2 iv.933350 v.430913 iv.705522 CREB3L1 three.622720 3.870238 3.514235 RPS11 ix.970795 ten.387306 10.124355 PNMA1 7.959072 8.356616 8.132725 > This data of ~12042 genes from Affymetrix HG-U133A is preprocessed and normalized using the RMA method, which results in profiles usually ranging from ii to 14 as shown in the histogram below (Fig.2):
Fig.2. Histogram of Affymetrix HG-U133A expression data.
Affymetrix Human Exon 1.0 ST Array
Affymetrix Man Exon 1.0 ST Array is another common platform to profile gene level expression. Unlike HG-U133A where probes were designed to embrace the entire genome, probes in Human Exon i.0 ST Array are designed to cover exonic regions of the genes.
# Get Affymetrix Human Exon 1.0 ST Array expression for ovarian cancer patients > huex.ov <- getTCGA(disease="OV", information.type="mRNA_Array", type="Huex") # Look at the information > dim(huex.ov$dat) [1] 18632 594 > caput(huex.ov$dat[,1:iii]) TCGA-01-0628-11A-01R-0361-03 TCGA-01-0630-11A-01R-0361-03 TCGA-01-0631-11A-01R-0361-03 C9orf152 5.684061 v.740597 5.379095 ELMO2 6.940721 6.955458 6.859219 RPS11 12.112295 11.603771 11.435075 CREB3L1 half dozen.141857 6.470617 six.821231 PNMA1 9.069748 9.193657 8.966602 MMP2 7.993179 8.551214 9.651369 This level-3 information is quantile-normalized, consisting of ~18632 genes with values range from 2 to 14 as shown in the post-obit histogram (Fig.3):
Fig.iii. Histogram of Affymetrix Human Exon one.0 ST Array information.
4.1.2. RNA-Seq gene expression information
TCGA also provides gene expression data measured via RNA-Sequencing applied science, particularly via Illumina Genome Analyzer (GA) or Illumina HiSeq 2000. For either platform, the data is available in two forms: RNA-Seq (preprocessed using the beginning pipeline), and RNA-SeqV2 (preprocessed using version two analysis).
RNA-Seq
Our package allows the user to download RNA-Seq information from both Illumina GA and Illumina HiSeq, depending on data availability. If gene expression is profiled via both platforms, merely the data from HiSeq will be imported. In addition, users have the following options for which type of data to import: raw-counts (type="count") or the normalized counts (blazon="RPKM").
# Get raw count, RNA-Seq data for ovarian cancer patients > counts.ov <- getTCGA(affliction="OV", data.blazon="RNASeq") # or equivalently > counts.ov <- getTCGA(illness="OV", data.type="RNASeq", blazon="count") # Look at the data: > dim(counts.ov$dat) [1] 19990 299 > head(counts.ov$dat[,one:3]) TCGA-04-1348-01A-01R-1565-13 TCGA-04-1357-01A-01R-1565-13 TCGA-04-1362-01A-01R-1565-xiii AADACL3 36 64 30 AADACL4 0 0 0 ABCA4 575 92 226 ABCB10 2308 1864 1572 ABCD3 4890 1437 5160 ABL2 3591 2187 3038 # Get normalized counts, from RNA-Seq for ovarian cancer patients > rnaseq.ov <- getTCGA(affliction="OV", data.type="RNASeq", blazon="RPKM") # Look at the information: > dim(rnaseq.ov$dat) [ane] 19990 299 > head(rnaseq.ov$dat[,ane:3]) TCGA-04-1348-01A-01R-1565-thirteen-2 TCGA-04-1357-01A-01R-1565-13-2 TCGA-04-1362-01A-01R-1565-13-ii AADACL3 0.0539 0.2341 0.0573 AADACL4 0.0000 0.0000 0.0000 ABCA4 0.4204 0.1647 0.2120 ABCB10 3.6155 seven.1614 3.1549 ABCD3 5.3981 3.8884 7.2949 ABL2 1.6614 ii.4809 1.7997 RNA-SeqV2
For RNA-Sequencing data from the 2d analysis pipeline (RNASeqV2), our package provides importation of data from Illumina HiSeq just equally this is the most mutual profiling method. The values returned are the RSEM values for the genes, which unremarkably contains 20501 genes, with RSEM values ranging from 0 to x^half-dozen.
# Go RSEM values, RNA-SeqV2 data for ovarian cancer patients rsem.ov <- getTCGA(disease="OV", information.type="RNASeq2") # Look at the data > dim(rsem.ov$dat) [1] 20501 307 > head(rsem.ov$dat[,i:3]) TCGA-04-1348-01A-01R-1565-thirteen TCGA-04-1357-01A-01R-1565-xiii TCGA-04-1362-01A-01R-1565-13 A1BG 66.4695 65.5664 41.6412 A1CF 0.0000 0.0000 0.3310 A2BP1 0.2689 0.6510 4.3025 A2LD1 221.5219 141.2826 265.8161 A2ML1 7.5289 54.6875 5.6263 A2M 5899.8279 9384.4401 3350.4207 4.two. Import miRNA expression
Similarly to gene expression, expression profiles of micro-RNAs (miRNA) are available from both the microarray platforms and the sequencing platforms. However, while the miRNA-Seq information is available for all cancer types, miRNA-array information is but available for three cancer types: GBM, GBMLGG, and OV.
iv.2.1. miRNA-assortment data
The miRNA expression data available for download with our packet is the level-iii Altitude Weighted Discrimination (DWD) batch adjusted miRNA expression profiled from Agilent 8 x 15K Human being miRNA-specific microarray (H-miRNA_8x15K) or Agilent Man miRNA Microarray Rel12.0 (H-miRNA_8x15Kv2 - for OV only).
# Get miRNA expression data via microarray for ovarian cancer patients mir.ov <- getTCGA(disease="OV", data.type="miRNA_Array") # Wait at the data: > str(mir.ov) Listing of 3 $ dat : num [i:799, one:594] 4.67 iv.58 iv.62 4.75 4.55 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:799] "ebv-miR-BART1-3p" "ebv-miR-BART1-5p" "ebv-miR-BART10" "ebv-miR-BART10*" ... .. ..$ : chr [one:594] "TCGA-01-0628-11A-01T-0364-07" "TCGA-01-0630-11A-01T-0364-07" "TCGA-01-0631-11A-01T-0364-07" "TCGA-01-0633-11A-01T-0364-07" ... $ clinical : NULL $ merged.dat: NULL > head(mir.ov$dat[,1:3]) TCGA-01-0628-11A-01T-0364-07 TCGA-01-0630-11A-01T-0364-07 TCGA-01-0631-11A-01T-0364-07 ebv-miR-BART1-3p 4.668551 4.640936 v.170681 ebv-miR-BART1-5p iv.579950 iv.696028 4.855175 ebv-miR-BART10 four.618877 iv.684807 4.742815 ebv-miR-BART10* 4.747920 4.730715 4.734904 ebv-miR-BART11-3p 4.553927 4.648465 4.605977 ebv-miR-BART11-5p 4.698162 iv.642823 4.711533 iv.2.ii. miRNA-Seq data
miRNA expression was profiled for well-nigh cancer types using Illumina HiSeq 2000 miRNA Sequencing. Similar to RNA-Seq, two types of miRNA-Seq data tin can be downloaded: raw counts (blazon="count") or RPMMM (type="rpmmm"). Examples below prove how to obtain each type of miRNA-Seq data:
# Go miRNA expression read counts for ovarian cancer patients mirseq.ov <- getTCGA(disease="OV", data.blazon="miRNASeq") # or equivalently mirseq.ov <- getTCGA(disease="OV", information.blazon="miRNASeq", type="count") # # Get miRNA expression RPMMM values for ovarian cancer patients mirseqn.ov <- getTCGA(affliction="OV", data.type="miRNASeq", type="rpmmm") > str(mirseq.ov) List of three $ dat : int [1:705, i:461] 124097 247488 124173 523763 133064 39686 31230 47 11523 1375 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:705] "hsa-allow-7a-1" "hsa-permit-7a-2" "hsa-let-7a-3" "hsa-let-7b" ... .. ..$ : chr [1:461] "TCGA-04-1331-01A-01R-1569-13" "TCGA-04-1332-01A-01R-1564-13" "TCGA-04-1336-01A-01R-1564-13" "TCGA-04-1337-01A-01R-1564-xiii" ... $ clinical : NULL $ merged.dat: NULL > head(mirseq.ov$dat[,1:3]) TCGA-04-1331-01A-01R-1569-thirteen TCGA-04-1332-01A-01R-1564-13 TCGA-04-1336-01A-01R-1564-13 hsa-let-7a-1 124097 9682 3037 hsa-permit-7a-2 247488 19634 6274 hsa-permit-7a-iii 124173 9825 3132 hsa-allow-7b 523763 88962 66321 hsa-allow-7c 133064 22926 5265 hsa-let-7d 39686 8955 14718 > caput(mirseqn.ov$dat[,1:3]) TCGA-04-1331-01A-01R-1569-13-one TCGA-04-1332-01A-01R-1564-13-ane TCGA-04-1336-01A-01R-1564-thirteen-1 hsa-let-7a-ane 30713.962 11269.88 3953.866 hsa-let-7a-ii 61253.190 22854.04 8168.112 hsa-let-7a-3 30732.772 11436.33 4077.546 hsa-let-7b 129631.151 103552.07 86343.214 hsa-let-7c 32933.291 26685.94 6854.496 hsa-let-7d 9822.271 10423.65 19161.343 4.3. Import mutation information
The level-three files of called mutations from the Wide Firehose are saved in Mutation Notation Format (MAF). Each MAF file lists mutations found for the particular patient. Since each patient may accept dissimilar mutations, the number of genes in each MAF file vary from patient to patient. Hence, the MAF files are not in a format ready for statistical analysis. To solve this, our package aggregates the obtained MAF files into a matrix of dimension gene x patient, with a value of one in cell(i,j) if a mutation is found in gene-i from patient-j, and a naught otherwise.
Our package permits two options for the mutation data:
- Only somatic non-silent mutations called (
blazon="somatic") which is the default. This will but return those mutations with "somatic"" equally the mutation-condition and not with "silent" as the variant-nomenclature in the original MAF files. - All mutations called, (
type="all").
Here are some example on downloading somatic mutation data:
# Get somatic not-silent mutations for ovarian cancer patients mut.ov <- getTCGA(disease="OV", data.blazon="Mutation") # or equivalently mut.ov <- getTCGA(affliction="OV", data.type="Mutation", type="somatic") > str(mut.ov) List of 3 $ dat : num [one:5821, 1:232] one one ane 1 1 ane 1 0 0 0 ... ..- attr(*, "dimnames")=Listing of two .. ..$ : chr [1:5821] "PPP2R5A" "TBC1D4" "ICAM1" "CPS1" ... .. ..$ : chr [1:232] "TCGA-24-1471" "TCGA-thirteen-0899" "TCGA-23-1021" "TCGA-04-1347" ... $ clinical : NULL $ merged.dat: Zilch > head(mut.ov$dat[,1:six]) TCGA-24-1471 TCGA-13-0899 TCGA-23-1021 TCGA-04-1347 TCGA-ten-0930 TCGA-23-1027 PPP2R5A one 0 0 0 0 0 TBC1D4 1 0 0 0 0 0 ICAM1 1 0 0 0 0 0 CPS1 1 0 0 0 0 0 EP300 1 0 0 0 0 0 PEBP4 ane 0 0 0 0 0
Note that in the example below, the number of genes is much larger if 1 is importing all mutations; only the information matrix is considerably sparser.
# Go all mutations chosen for ovarian cancer patients > allmut.ov <- getTCGA(illness="OV", information.type="Mutation", blazon="all") > str(allmut.ov) Listing of three $ dat : num [1:10060, 1:316] 1 1 ane 1 1 1 1 ane 1 1 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:10060] "HRNR" "LHX9" "PPP2R5A" "ZP4" ... .. ..$ : chr [ane:316] "TCGA-24-1471" "TCGA-xiii-0899" "TCGA-23-1021" "TCGA-04-1347" ... $ clinical : NULL $ merged.dat: Zip 4.4. Import methylation expression
Dna methylation profiles were obtained either via Illumina Infinium HumanMethylation27 BeadChip or the HumanMethylation450 BeadChip. The former platform probes for ~27,000 CpG sites (thus commonly known 27K) and the latter probes for ~450,000 CpG sites. Our parcel allows users to import either methylation information when available, via type="27K" (default) or type="450K". For case:
# Get 27K methylation profiles for ovarian cancer patients methyl.ov <- getTCGA(disease="OV", data.type="Methylation", type="27K") # or equivalently methyl.ov <- getTCGA(disease="OV", information.type="Methylation") # Get 450K methylation profiles for ovarian cancer patients methyl45.ov <- getTCGA(disease="OV", data.blazon="Methylation", type="450K") Note that the matrix of the methylation profiles returned is Non aggregated at the gene level. Each row in the data matrix (dat) is a probe from the methylation assay, which represents a CpG site. Since genes ofttimes contain more than one CpG site and each CpG site can differ significantly in the methylation level, gene level assemblage is less desirable. Hence, our bundle returns the probe-level data.
# Expect at the data > str(methyl.ov) List of 4 $ dat : num [1:27578, 1:612] 0.7994 0.339 0.0281 0.6012 NA ... ..- attr(*, "dimnames")=Listing of two .. ..$ : chr [1:27578] "cg00000292" "cg00002426" "cg00003994" "cg00005847" ... .. ..$ : chr [i:612] "TCGA-01-0628-11A-01D-0383-05" "TCGA-01-0630-11A-01D-0383-05" "TCGA-01-0631-11A-01D-0383-05" "TCGA-01-0633-11A-01D-0383-05" ... $ cpgs :'information.frame': 27578 obs. of 3 variables: ..$ Gene_Symbol : chr [i:27578] "ATP2A1" "SLMAP" "MEOX2" "HOXD3" ... ..$ Chromosome : chr [1:27578] "sixteen" "iii" "7" "two" ... ..$ Genomic_Coordinate: int [ane:27578] 28890100 57743543 15725862 177029073 148822837 93862594 93813777 11980953 89290921 0 ... $ clinical : NULL $ merged.dat: Naught > caput(methyl.ov$dat[,one:3]) TCGA-01-0628-11A-01D-0383-05 TCGA-01-0630-11A-01D-0383-05 TCGA-01-0631-11A-01D-0383-05 cg00000292 0.79940858 0.62039417 0.91148841 cg00002426 0.33900444 0.18030460 0.29385049 cg00003994 0.02811930 0.03607298 0.06159123 cg00005847 0.60116497 0.64955777 0.47394908 cg00006414 NA NA NA cg00007981 0.01881682 0.01803597 0.02904295 > head(methyl.ov$cpgs) Gene_Symbol Chromosome Genomic_Coordinate cg00000292 ATP2A1 16 28890100 cg00002426 SLMAP 3 57743543 cg00003994 MEOX2 7 15725862 cg00005847 HOXD3 2 177029073 cg00006414 ZNF425;ZNF398 vii 148822837 cg00007981 PANX1 11 93862594 As shown above, the returned listing for methylation data has one additional chemical element, cpgs, which contains the probe's genomic information. Each row of cpgs is a CpG site and contains cistron symbols and genomic coordinates in columns; rows in cpgs are of the same number and order every bit the methylation profiles data dat.
Values in the methylation information are the beta values. They range from zero to one (as shown in the histogram in Fig.iv below), with values greater than 0.v representing hypermethylation, values less than 0.5 representing hypomethylation.
Fig.iv. Histogram of methylation data.
iv.5. Import re-create number changes
2 major high-throughput platforms used to profile copy number changes are SNP-arrays and array comparative genomic hybridization (aCGH). Our parcel permits users to download both of these information types, although only SNP-array information is available for nigh cancer types. The copy number changes downloaded with our package are aggregated at the gene level; specifically, we merge the level-3 segments provided from TCGA (with reference genome build hg19) at the gene level. Values in the dat matrix of the returned listing object are the log-ratio of copy number between samples to reference samples. Thus, a positive log-ratio indicates a copy number gain and negative value indicates a DNA copy number loss.
Two forms of copy number changes profiled using SNP-arrays (Affymetrix Genome-Wide Human SNP Array 6.0) are bachelor: information.type="CNA_SNP" gives somatic copy number alterations (CNAs), which are somatic changes to chromosome structure that result in gains or losses in copies of sections of Dna; information.type="CNV_SNP" gives copy number variants (CNVs), which are the re-create number ratios from somatic cells minus the ratios from germlines. Examples of these are given below:
# Get SNP array CNA for ovarian cancer patients cnasnp <- getTCGA(disease="OV", information.type="CNA_SNP") # or # Go SNP array CNV for ovarian cancer patients cnvsnp <- getTCGA(disease="OV", data.type="CNV_SNP") # Await at the data > head(cnvsnp$dat[,1:3]) TCGA-01-0628-11A-01D-0356-01 TCGA-01-0630-11A-01D-0356-01 TCGA-01-0631-11A-01D-0356-01 A1BG -0.0043 0.0002 -0.0797 A2M 0.0028 -0.0003 0.k A2MP1 0.0028 -0.0003 0.thou NAT1 0.0047 0.0050 -0.0045 NAT2 0.0047 0.0050 -0.0045 SERPINA3 0.0024 0.0031 0.0501 Re-create number changes profiled via assortment-CGH are measured with 2 different types of arrays: Agilent Human being Genome CGH Microarray 244A (type="244A") and Agilent Human Genome CGH Custom Microarray 2x415K (blazon="415K"). Unlike CNA/CNV profiled with SNP-arrays that are bachelor for almost cancer types, aCGH-profiled CNA data are bachelor for but a few cancer types and most of them just for Custom Array 2x415. Examples of how to obtain this information are given below:
# Go aCGH CNA for ovarian cancer patients cnacgh <- getTCGA(disease="OV", data.blazon="CNA_CGH") # or equivalently cnacgh <- getTCGA(disease="OV", data.type="CNA_CGH", type="415K") # Expect at the information > names(cnacgh) [1] "dat" "clinical" "merged.dat" > caput(cnacgh$dat[,1:3]) TCGA-04-1353-01A-01D-1046-02 TCGA-04-1353-11B-01D-1046-02 TCGA-04-1369-01A-02D-1046-02 A1BG -0.7638 -0.2282 -0.5277 A2M -0.0284 -0.0774 0.0367 A2MP1 -0.0284 -0.0774 0.0367 NAT1 -0.8137 -0.1329 -one.1487 NAT2 -0.8137 -0.1329 -1.1487 SERPINA3 -0.2948 -0.0776 -0.1730 While merging the segments of copy number changes at cistron-level, nosotros filter out chromosome Y past default, which means the to a higher place three commands are equivalent to the following:
cnasnp <- getTCGA(disease="OV", information.blazon="CNA_SNP", filter="Y") cnvsnp <- getTCGA(disease="OV", data.type="CNV_SNP", filter="Y") cnacgh <- getTCGA(illness="OV", data.blazon="CNA_CGH", type="415K", filter="Y") To include all genes, give an empty string to filter, for example:
cnacgh.gbm <- getTCGA(disease="GBM", data.type="CNA_CGH", type="415K", filter="") Since in that location are certain factor symbols shared by two chromosomal locations in the sex activity chromosomes, we suffixed these genes past the chromosome to eliminate ambiguities. For example, in the cnacgh.gbm$dat object obtained from the above command, "SLC25A6_X" and "SLC25A6_Y" stand for SLC25A6 genes in chromosome 10 and SLC25A6 in chromosome Y respectively.
Similarly, by providing a vector of characters with more than one chromosome to filter, the function will filter out any genes belonging to the specified chromosomes. For example, to get CNA autosomal chromosomes simply (excluding sex chromosomes):
cnacgh <- getTCGA(disease="OV", data.blazon="CNA_CGH", filter=c("Ten", "Y")) 5. Miscellaneous Functions
5.one. Combine multiple OMICs profiles
Our package includes a role named OMICSBind that allows users to combine ii matrices of omics-profiles for the same set of patients into a single object. In order to combine more than two types of omics-profiles together, users tin can phone call the OMICSBind office multiple times. In the case beneath, we testify how to combine the RNA-Seq data from the 2d pipeline analysis, methylation profiles from 27K, and mutation information for ovarian cancer patients:
# Become the RNA-Seq, methylation, and mutation profiles for OV cancer patients seq <- getTCGA(affliction="OV", information.type="RNASeq2") meth <- getTCGA(disease="OV", data.type="Methylation", type="27K") mut <- getTCGA(disease="OV", data.type="Mutation", blazon="all") # Now, merge the three OMICs-data into ane R object # step ane: merge RNA-Seq and mutation data m1 <- OMICSBind(dat1 = seq$dat, dat2 = mut$dat) # stride ii: further concatenate the methylation data to the merged data-object m2 <- OMICSBind(dat1 = m1$merged.information, dat2 = meth$dat) # Look at the data > dim(seq$dat) [i] 20501 307 > dim(mut$dat) [1] 10060 316 > dim(meth$dat) [1] 27578 612 > dim(m1$merged.data) [1] 30561 185 > dim(m2$merged.data) [i] 58139 185 Every bit shown above, the final merged information matrix m2$merged.information has 185 samples. These are the tumor samples of patients common in RNA-Seq, mutation, and methylation information. On the other hand, m2$merged.data has 58139 rows, which is the combination of 20501 genes from RNA-Seq, 10060 genes from mutation, and 30561 CpG sites from methylation information. In order to distinguish factor symbols from different omics-profile, we affix the cistron names with d1 and d2 in the merged data. Specifically, d1 is affixed to gene names of the beginning omics-information input to OMICSBind function, d2 is affixed to the names of the second input object:
> str(m1) List of iii $ merged.information: num [1:30561, 1:185] 66.469 0 0.269 221.522 7.529 ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [ane:30561] "d1.A1BG" "d1.A1CF" "d1.A2BP1" "d1.A2LD1" ... .. ..$ : chr [1:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... $ X : num [1:20501, 1:185] 66.469 0 0.269 221.522 7.529 ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [1:20501] "d1.A1BG" "d1.A1CF" "d1.A2BP1" "d1.A2LD1" ... .. ..$ : chr [1:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... $ Y : num [ane:10060, ane:185] 0 0 0 0 0 0 0 0 0 0 ... ..- attr(*, "dimnames")=Listing of two .. ..$ : chr [1:10060] "d2.HRNR" "d2.LHX9" "d2.PPP2R5A" "d2.ZP4" ... .. ..$ : chr [1:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... > str(m2) Listing of 3 $ merged.data: num [1:58139, 1:185] 66.469 0 0.269 221.522 vii.529 ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [one:58139] "d1.d1.A1BG" "d1.d1.A1CF" "d1.d1.A2BP1" "d1.d1.A2LD1" ... .. ..$ : chr [one:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... $ Ten : num [i:30561, 1:185] 66.469 0 0.269 221.522 7.529 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:30561] "d1.d1.A1BG" "d1.d1.A1CF" "d1.d1.A2BP1" "d1.d1.A2LD1" ... .. ..$ : chr [1:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... $ Y : num [1:27578, one:185] 0.7991 0.0968 0.0218 0.8795 NA ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [ane:27578] "d2.cg00000292" "d2.cg00002426" "d2.cg00003994" "d2.cg00005847" ... .. ..$ : chr [1:185] "TCGA-04-1348" "TCGA-04-1357" "TCGA-04-1362" "TCGA-04-1364" ... five.two. Separating OMICs profiles past sample types
The TCGA OMICs profiles downloaded for a cancer type often include samples of different types, such equally from tumor tissue, normal tissue, or other controls. The type of sample is encoded in the sample's BCR code. Users who are interested in studying just tumor samples, for example, would need to parse the codes to excerpt only these samples. To spare users the difficulties of decoding the sample ID (BCR code) to extract particular types of samples, our package includes a role, SampleSplit to accomplish this job in a user-friendly manner. Specifically, SampleSplit volition take the omics data matrix equally input and split up the matrix according to the sample types. The object returned is a list of iii matrices corresponding to the omics profiles of principal solid tumor, recurrent solid tumor, and normal tissues/blood samples.
# Get the RNA-Seq V2 profiles for LUSC cancer patients lusc.rnaseq2 <- getTCGA(disease="LUSC", data.type="RNASeq2") # Split the OMICs data past sample types lusc.rnaseq2.bytype <- SampleSplit(lusc.rnaseq2$dat) # Expect at the data > str(lusc.rnaseq2.bytype) Listing of 3 $ primary.tumor : num [1:20501, i:501] 742 0 0 170 128 ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [1:20501] "A1BG" "A1CF" "A2BP1" "A2LD1" ... .. ..$ : chr [1:501] "TCGA-18-3406-01A-01R-0980-07" "TCGA-18-3407-01A-01R-0980-07" "TCGA-18-3408-01A-01R-0980-07" "TCGA-18-3409-01A-01R-0980-07" ... $ recurrent.tumor: num[i:20501, 0 ] ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [ane:20501] "A1BG" "A1CF" "A2BP1" "A2LD1" ... .. ..$ : NULL $ normal : num [one:20501, one:51] 49.49 0 1.04 105.22 0 ... ..- attr(*, "dimnames")=Listing of 2 .. ..$ : chr [1:20501] "A1BG" "A1CF" "A2BP1" "A2LD1" ... .. ..$ : chr [1:51] "TCGA-22-4593-11A-01R-1820-07" "TCGA-22-4609-11A-01R-2125-07" "TCGA-22-5471-11A-01R-1635-07" "TCGA-22-5472-11A-11R-1635-07" ... 5.3. Getting matched tumor-normal OMICs profiles
Our package also includes a utility function, TumorNormalMatch to prepare the downloaded omics-profiles for pairwise assay between tumor tissues and normal tissues/blood samples in a user-friendly manner. In the example below, we show how to starting time go the RNA-Seq data of assay pipeline ii from LUSC patients and so split the data into tumor tissues and normal tissues/blood samples for the same set of patients.
# Get RNA-SeqV2 information for LUSC patients > lusc.rnaseq2 <- getTCGA(affliction="LUSC", data.blazon="RNASeq2") # tumor-normal matched profiles > lusc.rnaseq2.tum.norm <- TumorNormalMatch(lusc.rnaseq2$dat) This command will render a list of two omics-profile data matrices with dimension of gene x patient. Both matrices accept the same dimension and order for both columns (patient) and rows (factor). Thus, enabling users to make direct pairwise comparisons.
# Wait at the information > str(lusc.rnaseq2.tum.norm) List of 2 $ primary.tumor: num [1:20501, 1:51] 94.67 0 iii.46 38.73 209.32 ... ..- attr(*, "dimnames")=Listing of two .. ..$ : chr [1:20501] "A1BG" "A1CF" "A2BP1" "A2LD1" ... .. ..$ : chr [one:51] "TCGA-22-4593" "TCGA-22-4609" "TCGA-22-5471" "TCGA-22-5472" ... $ normal : num [1:20501, 1:51] 49.49 0 1.04 105.22 0 ... ..- attr(*, "dimnames")=Listing of two .. ..$ : chr [ane:20501] "A1BG" "A1CF" "A2BP1" "A2LD1" ... .. ..$ : chr [i:51] "TCGA-22-4593" "TCGA-22-4609" "TCGA-22-5471" "TCGA-22-5472" ... 6. Example: Approved correlation analysis
The major objective of developing TCGA2STAT packet is to provide statisticians with a gear up of user friendly tools to obtain and prepare omics-profiles from TCGA cancer patients for statistical analysis. In this example, we show how TCGA2STAT can enable easy and convenient integrated assay of 2 types of TCGA omics profiles using canonical correlation analysis.
First, we import methylation profiles and RNA-Seq cistron expression data for LUSC patients via the getTCGA office from our package:
# Get data > lusc.methyl <- getTCGA(illness="LUSC", data.type="Methylation") > lusc.rnaseq2 <- getTCGA(disease="LUSC", data.type="RNASeq2", clinical=Truthful) To reduce computation fourth dimension, nosotros only include CpG sites and genes that are at the pinnacle 1% of variance across patients:
# Filter data > met.var <- employ(lusc.methyl$dat, 1, var) > met.data <- subset(lusc.methyl$dat, met.var >= quantile(met.var, 0.99, na.rm=T) & !is.na(met.var)) > rnaseq2.var <- employ(log10(1+lusc.rnaseq2$dat), ane, var) > rnaseq.data <- subset(log10(one+lusc.rnaseq2$dat), rnaseq2.var >= quantile(rnaseq2.var, 0.99, na.rm=T) & !is.na(rnaseq2.var)) This will give us 234 CpG sites and 206 genes. Since canonical correlation analysis requires two unlike types of features measured on the aforementioned set of samples, we employ the OMICSBind part in the parcel to get R objects with methylation profiles and gene expression profiles for the same set of LUSC patients:
> met.rnaseq2 <- OMICSBind(dat1 = rnaseq.data, dat2= met.data) # Look at the data > str(met.rnaseq2) List of 3 $ merged.data: num [1:440, 1:130] ii.112 1.654 0.252 2.874 1.391 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:440] "d1.A2ML1" "d1.ADAM23" "d1.ADCY8" "d1.ADH1B" ... .. ..$ : chr [1:130] "TCGA-18-3406" "TCGA-18-3407" "TCGA-18-3408" "TCGA-eighteen-3409" ... $ X : num [1:206, one:130] 2.112 1.654 0.252 ii.874 1.391 ... ..- attr(*, "dimnames")=Listing of ii .. ..$ : chr [i:206] "d1.A2ML1" "d1.ADAM23" "d1.ADCY8" "d1.ADH1B" ... .. ..$ : chr [1:130] "TCGA-18-3406" "TCGA-18-3407" "TCGA-18-3408" "TCGA-18-3409" ... $ Y : num [1:234, 1:130] 0.7561 0.6541 0.642 0.024 0.0258 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [i:234] "d2.cg00061629" "d2.cg00202711" "d2.cg00332153" "d2.cg00463577" ... .. ..$ : chr [1:130] "TCGA-18-3406" "TCGA-18-3407" "TCGA-18-3408" "TCGA-18-3409" ... Now that the data is ready, we acquit out canonical correlation analysis using the CCA package. Nosotros utilize regularized approved correlation analysis as the number of features for both sets of omics-profiles are larger than the number of samples.
library(CCA) # run regularized-cca > lusc.cc <- rcc(t(met.rnaseq2$X), t(met.rnaseq2$Y), 0.75025, 0.5005) # compute the canonical loadings > lusc.cc2 <- comput(t(met.rnaseq2$X), t(met.rnaseq2$Y), lusc.cc) We tin can visualize the results of the assay. In order to generate plots shown in our paper, nosotros wrote plotting functions that can be institute here: script. Showtime, we visualize how the factor expression profiles correlate with the methylation profiles:
> plotNice2(lusc.cc2, d1=1, d2=ii, XY="Ten", cex=0.7) 
Fig.5. Canonical loadings of genes.
(Note - fill in Xxx) From the plot in Figure five, we can see that gene XIST is highly correlated with the set up of CpG sites. Similarly, we can run into how methylation profiles correlate with the gene expression profiles:
> plotNice2(lusc.cc2, d1=i, d2=two, XY="Y", cex=0.7) 
Fig.6. Approved loadings of CpG sites.
With approved correlation analysis, we can also investigate patterns among patients by projecting onto the approved loadings. Here, we are interested in determining whether the estimated canonical loadings separate patients into distinct groups or cancer subtypes. Start, nosotros project onto the offset ii canonical loadings and visualize the patients below:
> plt.indiv(lusc.cc, d1=1, d2=2, ind.names=rep("+", nrow(t(met.rnaseq2$Ten)))) 
Fig.7. Patient clusters.
In Effigy 7, nosotros can see that patients are well-clustered into two groups. We re-plot this with two different colors in order to clearly see the two groups:
score <- lusc.cc2$xscores[,1:2] # Parameters were estimated based on visual inspection grp1 <- rownames(score)[(score[,1]<0 & score[,two]<0) | (score[,1]>0 & score[,2]< -1.5)] grp2 <- setdiff(rownames(score), grp1) grp <- ifelse(rownames(score) %in% grp1, i, ii) colors = colors()[c(153,220)] plotNice2.indiv(lusc.cc2, d1=ane, d2=ii, cols=colors[grp], cex=ii.five,pchs=c(twenty,20)[grp]) 
Fig.8. Nicer plot of patient clusters as in Fig.7.
Next, we investigate whether the two patient groups found by our assay are clinically meaningful past seeing whether the groups exhibit differing survival outcomes. To this end, we perform a Kaplan-Meier analysis on the overall patient survival time to determine whether the survival probabilities betwixt the two identified patient subgroups are significantly different:
# Get survival data from the imported RNASeq2 object os.score <- lusc.rnaseq2$merged.dat[,1:iii] # prepare the information to make sure they match the ii clusters from above os <- os.score[,-1] rownames(os) <- bone.score[,one] os <- os[rownames(lusc.cc2$xscores),] plotdat <- information.frame(grp, os) plotdat <- plotdat[order(plotdat[,1]),] # KM-plot mypamr.plotsurvival(plotdat[,i], plotdat[,3]/365, plotdat[,2], c(colors()[c(153,220)]), lty=c(2,1), lwd=c(2,two)) 
Fig.nine. KM-plot of patients' overall survival.
From the survival plot in Figure ix, we tin can conclude that the canonical loadings of gene expression and methylation profiles divide patients into two clinically distinct subgroups: the 'black group' is the low-adventure group and the 'greyness group' is the high-risk group. The divergence in overall survival between these ii groups is border-line statistically pregnant (log-rank p-value < 0.06).
7. Appendix
vii.ane. Appendix A
Summary of cancer types and omics-profiles. Go File
7.ii. Appendix B
Values permitted for data.blazon and type in getTCGA function. Get File
7.3. Appendix C
Summary of clinical variables common in near cancer types. Go File
DOWNLOAD HERE
Posted by: ringfavy1940.blogspot.com
