[GRanges] How to append query metadata to subjectHits
Hi all,
I'm I found overlaps between couple of genes and TCGA copy number calls but I've been unsuccessful in appendding metadata (hgnc_symbol) to subjectHits
data look like this; 1. sample of genes of interest
> my.cordinates.hg19.gr
GRanges object with 3 ranges and 3 metadata columns:
seqnames ranges strand | hgnc_symbol band gene_biotype
<Rle> <IRanges> <Rle> | <character> <character> <character>
[1] 8 128747680-128753674 * | MYC q24.21 protein_coding
[2] 5 1253262-1295184 * | TERT p15.33 protein_coding
[3] 17 7565097-7590856 * | TP53 p13.1 protein_coding
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
TCGA absolute copy calls in Granges;
> Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities
GRanges object with 159857 ranges and 21 metadata columns:
seqnames ranges strand | Sample Num_Probes Length Modal_HSCN_1 Modal_HSCN_2 Modal_Total_CN Subclonal_HSCN_a1 Subclonal_HSCN_a2 Cancer_cell_frac_a1 Ccf_ci95_low_a1
<Rle> <IRanges> <Rle> | <character> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
[1] 1 564621-1510801 * | TCGA-3C-AAAU_WHITE 53 946180 1 2 3 0 1 0.25 0.08383
[2] 1 1688192-16142960 * | TCGA-3C-AAAU_WHITE 4180 14454768 1 2 3 0 1 0.06 0.02732
[3] 1 16165661-63472103 * | TCGA-3C-AAAU_WHITE 13591 47306442 1 2 3 0 1 0.04 0.00615
[4] 1 63472868-72759524 * | TCGA-3C-AAAU_WHITE 3196 9286656 1 2 3 0 0 0.03 0.00000
[5] 1 72811904-85632596 * | TCGA-3C-AAAU_WHITE 4069 12820692 1 2 3 0 0 0.04 0.00292
... ... ... ... . ... ... ... ... ... ... ... ... ... ...
[159853] 21 28284150-39907857 * | TCGA-Z7-A8R6_WHITE 4229 11623707 1 1 2 0 0 0.00 0.00000
[159854] 21 39908195-48084820 * | TCGA-Z7-A8R6_WHITE 3214 8176625 1 1 2 0 0 0.05 0.01384
[159855] 22 16055207-24329711 * | TCGA-Z7-A8R6_WHITE 1878 8274504 1 1 2 0 0 0.05 0.00946
[159856] 22 24402321-44489022 * | TCGA-Z7-A8R6_WHITE 6479 20086701 1 1 2 0 0 0.04 0.00644
[159857] 22 44493823-51199978 * | TCGA-Z7-A8R6_WHITE 2807 6706155 1 1 2 0 0 0.06 0.02714
Ccf_ci95_high_a1 Cancer_cell_frac_a2 Ccf_ci95_low_a2 Ccf_ci95_high_a2 LOH Homozygous_deletion solution Type race ethnicity
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <character> <character> <character> <character>
[1] 0.39979 0.74 0.58132 0.89771 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[2] 0.08989 0.85 0.81351 0.87773 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[3] 0.06643 0.91 0.87122 0.92977 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[4] 0.05842 0.03 0.00000 0.05981 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[5] 0.06738 0.06 0.02157 0.08639 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
... ... ... ... ... ... ... ... ... ... ...
[159853] 0.02834 0.00 0.00000 0.02834 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[159854] 0.06983 0.01 0.00000 0.03265 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[159855] 0.07308 0.05 0.00946 0.07308 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[159856] 0.05948 0.04 0.00644 0.05948 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
[159857] 0.08524 0.06 0.02714 0.08524 0 0 new BRCA WHITE NOT HISPANIC OR LATINO
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
>
find overlaps between TCGA asbolute copy calls and ancestry- related genes
overlaps.tcga.ancestry.genes <- findOverlaps(query=my.cordinates.hg19.gr, subject=Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities, type = "within") ##complete overlaps
create empty column to append gene names for easy filtering for downstream analysis
Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities$Ancestry_genes <- NA
> str(Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities$Ancestry_genes) ##sanity checks
logi [1:159857] NA NA NA NA NA NA ...
extract indexes of subject hits
sub_hits.overlaps <- as.data.frame(Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities[subjectHits(overlaps.tcga.ancestry.genes)])
str(sub_hits.overlaps$Ancestry_genes) ##sanity checks
append genes names to corresponding overlap indexes in subject hits
Granges.merged.df.tcga.mastercalls.abs.BRCA.ID.ethnicities[queryHits(overlaps.tcga.ancestry.genes)]$Ancestry_genes <- sub_hits.overlaps[,grepl("hgnc_symbol",names(sub_hits.overlaps))]
Error in do.call(`[<-`, args) : replacement has length zero
Thanks for your help
1 answer
I think it is better if you use the metadata columns in the object, with the mcols function. This function lets you access the DataFrame containing the metadata inside of the GRanges. For example, these 2 GRanges:
> library(GenomicRanges)
> gr <- GRanges(
+ seqnames = c("chr3", "chr4", "chr2"),
+ ranges = IRanges(start = c(101,103,102), end = c(111,113,112), names = head(letters, 3)),
+ hgnc_symbol = c("MYC","TERT","TP53"))
> gr
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | hgnc_symbol
<Rle> <IRanges> <Rle> | <character>
a chr3 101-111 * | MYC
b chr4 103-113 * | TERT
c chr2 102-112 * | TP53
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
>
> tcga.calls <- GRanges(
+ seqnames = c("chr2", "chr4", "chr3"),
+ ranges = IRanges(start = c(102,70,500), end = c(105,105,510), names = head(letters, 3)),
+ metadata = c("a","b","c"))
> tcga.calls
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | metadata
<Rle> <IRanges> <Rle> | <character>
a chr2 102-105 * | a
b chr4 70-105 * | b
c chr3 500-510 * | c
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
You find the overlaps, and you can check their positions:
> overlaps <- findOverlaps(query = gr, subject = tcga.calls)
> queryHits(overlaps)
[1] 2 3
> subjectHits(overlaps)
[1] 2 1
Now create your new metadata column in the subject:
> mcols(tcga.calls)$Ancestry_genes <- NA
> tcga.calls
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | metadata Ancestry_genes
<Rle> <IRanges> <Rle> | <character> <logical>
a chr2 102-105 * | a <NA>
b chr4 70-105 * | b <NA>
c chr3 500-510 * | c <NA>
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
And finally assign the matched values:
> mcols(tcga.calls)$Ancestry_genes[subjectHits(overlaps)] <- mcols(gr)$hgnc_symbol[queryHits(overlaps)]
>
> tcga.calls
GRanges object with 3 ranges and 2 metadata columns:
seqnames ranges strand | metadata Ancestry_genes
<Rle> <IRanges> <Rle> | <character> <character>
a chr2 102-105 * | a TP53
b chr4 70-105 * | b TERT
c chr3 500-510 * | c <NA>
-------
seqinfo: 3 sequences from an unspecified genome; no seqlengths
Log in to answer this question.
More posts like this
-
count number of GC for a given genomic ranges
written by Alewa
-
rtracklayer import.bed adds one to start range
written by Katarina •
-
Aggregate differentially expressed genes across all contrast results using DEVis
written by Ana •
-
Help needed
written by Nai •
-
single cell sequencing understanding the structure of a data
written by Mo
-
How to analyse TCGA CNV data with GAIA R package
written by gina.abdelaal •
-
Calculate pairwise Fst from Plink (.bim .bed .fam) files in R?
written by bntjoanne •
-
Filtering GFF annotations file
written by alslonik
-
How to find key apoptotic markers across subtypes of a colorectal cancer dataset?
written by bio94 •
-
Boxplot each row in a dataset in R?
written by bio94 •
