Question: closest-features equivalent in Bioconductor for finding closest gene(s)
0
gravatar for Dimitris Polychronopoulos
3.3 years ago by
United Kingdom

Dear all,

I would like to get from a GRanges object, where I have stored genomic coordinates of my differentially expressed genes (using biomart and mm9), the closest gene(s) for every DE gene in my list. Using closest-features in unix I could do it easily:

closest-features total_genes_gr.bed all_mouse_genes_gr.bed

and then merge the dataframes.

But I would like to be consistent and do everything in R. Ideally, I would like to have those closest features / genes as a metadata column. Could you please suggest any ways of doing this?

Best wishes,

Dimitris

 

 

ADD COMMENTlink modified 3.3 years ago by Michael Lawrence11k • written 3.3 years ago by Dimitris Polychronopoulos50
Answer: closest-features equivalent in Bioconductor for finding closest gene(s)
2
gravatar for Michael Lawrence
3.3 years ago by
United States
Michael Lawrence11k wrote:

Should be straight-forward with ?nearest.

ADD COMMENTlink written 3.3 years ago by Michael Lawrence11k
Answer: closest-features equivalent in Bioconductor for finding closest gene(s)
1
gravatar for Dimitris Polychronopoulos
3.3 years ago by
United Kingdom

Nearest will give me the actual differentially expressed gene each time if i put as query the granges object of my diff expressed genes and as subject all mouse genes. Unless there is a special argument that gives you the second closest hit i will end up with the same list of genes

ADD COMMENTlink written 3.3 years ago by Dimitris Polychronopoulos50

Could just call nearest with one argument, all the genes. It will ignore self-hits. After annotating, subset to the DE genes.

ADD REPLYlink written 3.3 years ago by Michael Lawrence11k

Many thanks Michael for your suggestion. I could not make it work,though, but I think i found a way using follow. More specifically:

genes_index <- follow(total_DE_genes_gr,all_mouse_genes_gr,select=c("all")) #This gives me the second "best hit"

nearest_genes <- all_mouse_genes_gr[subjectHits(genes_index)]  #This is the list of mouse genes close to my DE genes

When I am trying to create a metadata column, though, having those genes side by side, I get the following error:

total_DE_genes_gr$nearest_gene <- all_mouse_genes_gr[subjectHits(genes_index)]
Error in `[[<-`(`*tmp*`, name, value = <S4 object of class "GRanges">) :
  20945 elements in value to replace 20956 elements

Probably because not all of my DE genes have a nearest gene, and therefore NAs are introduced? Any ways to circumvent this?

Dimitris

 

 

ADD REPLYlink written 3.3 years ago by Dimitris Polychronopoulos50
Answer: closest-features equivalent in Bioconductor for finding closest gene(s)
0
gravatar for Michael Lawrence
3.3 years ago by
United States
Michael Lawrence11k wrote:

Find nearest gene to each gene:

total_DE_genes_gr$nearest_gene <- all_mouse_genes_gr[nearest(all_mouse_genes_gr)[match(total_DE_genes_gr, all_mouse_genes_gr)]]
ADD COMMENTlink written 3.3 years ago by Michael Lawrence11k

Thanks again, Michael. Maybe I miss something, but when I run this, it returns:

total_DE_genes_gr$nearest_gene <- all_mouse_genes_gr[nearest(all_mouse_genes_gr)[match(total_DE_genes_gr, all_mouse_genes_gr)]]
Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) :
  subscript contains NAs or out-of-bounds indices

 

 

 

ADD REPLYlink written 3.3 years ago by Dimitris Polychronopoulos50
1

Two possibilities:

1) You have genes in the DE set that are not in the total set.

2) There are genes all by themselves on a chromosome/strand.

Btw, you should pass ignore.strand=TRUE to nearest() unless you only want the nearest gene on the same strand.

There is no way to encode an NA GRanges element, so if there are cases where a gene has no neighbor, you will need to think about what to do.

ADD REPLYlink written 3.3 years ago by Michael Lawrence11k

Thanks Michael for the suggestions, it worked!

ADD REPLYlink written 3.3 years ago by Dimitris Polychronopoulos50
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 167 users visited in the last hour