r query_tblastn.R

Posted 2021-05-24

library(dplyr)
library(readr)
library(tidyr)
library(stringr)
library(reutils)
library(XML)
library(reutils)


tblastn <- function(fasta, db = "refseq_genomic") 
{
  query_name <- paste0(fasta, ".blast.txt")
  comm <- paste("export BLASTDB=/usr/local/share/blast; tblastn -query ", fasta,
                "-db", db, 
                "-task tblastn",
                "-max_target_seqs 20000",
                "-outfmt '6  sscinames scomnames staxids qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore'",
                "-remote",
                "-entrez_query txid6237[ORGN]",
                ">",
                query_name)
  print(comm)
  system(comm)
  
  r <- read_tsv( query_name, col_names = c("Species",
                                           "Name",
                                           "TaxID",
                                           "QueryID",
                                           "SubjectID",
                                           "Percent_Identity",
                                           "Alignment_Length",
                                           "Mismatches",
                                           "Gap_Openings",
                                           "Q.Start",
                                           "Q.End",
                                           "S.Start",
                                           "S.End",
                                           "E",
                                           "Bits") ) %>%
    separate(SubjectID, into = c("name_drop", "gi", "ref_drop","accession"), sep = "\\|", extra = "drop", convert = T) %>%
    dplyr::select(-name_drop, -ref_drop) %>%
    dplyr::mutate(Name = sapply(unlist(efetch(accession, db="nuccore", "docsum")['//Item[@Name="Title"]/text()']), xmlValue) ) %>%
    dplyr::rename(POS_Start = S.Start, POS_End = S.End) %>%
    dplyr::mutate(CHROM = str_match(Name, "chromosome ([A-Za-z0-9])")[,2]) %>%
    dplyr::select(Species, Name, TaxID, QueryID, CHROM, POS_Start, POS_End, accession,  everything()) 
}

results <- tblastn("~/Desktop/Sap.fasta")