blast filtering

[james.r.m....@gmail.com]

Hello, 

I was wondering whether if there is a way of filtering a gl object based on gl.blast. 

I want to be 100% sure that the SNPs I have align to a reference genome (samples are from decades old museum material).

I've blasted the original gl object to my reference genome and have a df of SNPs that have been filtered with blast. How do I now filter my original gl SNP dartR object by this blast df?

Thanks!

James 

[Jose Luis Mijangos]

Hi James,

One way to do it is by filtering on evalue or bitscore, see for example, the code below. You can get more information about these two values here:

https://www.metagenomics.wiki/tools/blast/evalue 

library(dartR)

res <- gl.blast(x= platypus.gl,
                ref_genome = "chr_20_mOrnAna1.pri.v4.fa")
               
res_2 <- gl.filter.locmetric(res,
                             metric = "evalue",
                             lower = min(res$other$loc.metrics$evalue,
                                         na.rm = TRUE),
                             upper = 1e-20)

Cheers,

Luis 

[james.r.m....@gmail.com]

Hi Luis, 

Many thanks for this. Unfortunately,  when I run gl.blast it the object isn't a genlight object, just a data frame. How can I get the results of gl.blast back into a genlight object for me to to filter. See attached screenshot below:

Thanks, 

James 

[kw k]

Hi all,

I am super glad to see this question! I have the same problem, and Luis's suggestion would not work because we would like to map the SNPs using gl.blast to particular reference genomes. The filtering cannot be done by simply filtering the e-value or bit score in the locmetric of the genelight object because they are NOT the same reference genome I provided to DArT and mapped by them.

After blasting, if we want to filter, the filtering by loci function gl.keep.loci() or gl.drop.loci() only works with the specification of loci names in the loc.list argument. However, the command only accepts the names obtained using the locName(), instead of CloneID/TrimmedSequence you can see on the SNPs dataset in CSV from DArT. So, James and I will need the loci name workable with gl.keep.loci(), of the mapped sequences from gl.blast. We will need to match the blasted query sequence back to their loci name.

Does anyone have insight on how to do it? Also, as gl.blast output a dataframe with 'qseq', is it equivalent to the 'TrimmedSequence' in the SNPs genelight object?

Thank you so much!

Cheers, Christy

[james.r.m....@gmail.com]

Hi Christy, 

Thanks for the detailed layout of our shared problem! A thought I had would be to try using dplyr to merge 'qseq' with sequences in the genlight object.However, for that the query seq (qseq) will need to match the TrimmedSequence in the SNP genlight object. However, for some queries gaps are introduced as part of the BLAST step which means that a 1 for 1 dplyr merge won't work. Open to other suggestions!

Cheers, 

James 

[Jose Luis Mijangos]

Hi James,

There was a bug in the function when using a dartR object.

I have fixed it. To use the updated version of the function, you should install the developing version of dartR as shown below:

library(devtools)
install_github("green-striped-gecko/dartR@dev")

After installing it, please try to run the function again. 

Cheers,

Luis 

[james.r.m....@gmail.com]

Hi Luis, 

That's fixed it! Thank you so much!!

Cheers

James

[james.r.m....@gmail.com]

Sorry to find a new thing! It appears that the outputs of the blast aren't being added to $other$loc.metrics. Its outputting as a dart object, but can't find the results of the blast alignment

Cheers, 

James

[kw k]

Hi James and Luis,

Thanks for fixing the bug, Finally it can produce a genelight output! BUT I have the same problem with James - the loc.metrics of the output genelight object does not contain the blasted metric (e.g. qseq, pident, evalue etc) at all. However, it worked perfectly when I tried with Luis' code.....

Luis, could you please help to fix the code again? So grateful for your help! :)

Cheers,

Christy

[Jose Luis Mijangos]

It is working now!!

Please install the developing version of dartR again and try the function:

library(devtools)
install_github("green-striped-gecko/dartR@dev")

Cheers,

Luis 

[kw k]

Hi Luis,

Lovely! Thank you so much, it worked!

Best,

Christy