5
$\begingroup$

I'm attempting to visualize the results of my BLAST search in a way similar to the graphical display of the distribution of blast hits from the web BLAST.

For example, from my BLAST search:

enter image description here

except I would like the graph to contain information including information for the distribution of all 6139 hits across my query sequence.

Since this is a large number of sequences to visualize a distribution, I am considering using a score of the number of hits on a certain region, and having an output more similar to:

enter image description here

If I were to run my blast using the BLAST+ command line tool, and read the results into R, which parts of the output should I plot in order to recreate the graphical display of the distribution? Ideally I would like to make a reusable object that can make an equivalent plot for any BLAST input.

Thanks.

$\endgroup$
5
  • $\begingroup$ How do you get the score of number of hits on regions? I think you should get the score first and plot the score and region in R. If you want aliginment visulization, maybe this link help $\endgroup$ – l0o0 Jul 27 '17 at 3:53
  • 1
    $\begingroup$ What format would you use for the blast output? One of the standard ones, or are you open to defining your own? And would this be plotting the number of hits at each position of your query or their score? Or maybe their overall identity? What exactly is "Conservation" on your example axis? $\endgroup$ – terdon Jul 27 '17 at 9:06
  • $\begingroup$ @terdon I'm not wedded to any specific type of blast output, and am open to defining my own for this problem. I think that I would choose to plot the number of hits at each position instead of score. this link provides a description of the method used to make the conservation plot in my post in detail. $\endgroup$ – bluescholar1212 Jul 27 '17 at 16:32
  • $\begingroup$ But the number of hits at a given position is useless information if you don't take into account the quality of the relevant HSPs. You seem to want to plot conservation, if so, that's the wrong way to do it. For instance, you could have hundreds of hits at position 12, but if those hundreds are all tiny and very badly conserved, they are irrelevant. $\endgroup$ – terdon Jul 27 '17 at 16:35
  • $\begingroup$ I don't want to plot conservation. I want to know which regions of a given human protein are the most similar to bacterial sequences. If that is my question, does size/conservation necessarily matter? $\endgroup$ – bluescholar1212 Jul 27 '17 at 16:41
4
$\begingroup$

You probably want to include query start (qstart) and query end (qend) in your blast output.

Something like this:

blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'

In R you can take the "qstart:qend" from each line for density plot.

There are many ways in R to plot the densities of these start and end amino acids.

Let me show an example with a small data frame:

qstart <- c(200, 300, 250, 400, 500)
qend <- c(300, 450, 400, 600, 650)

df <- as.data.frame(cbind(qstart,qend))

aa <- vector()
i=1
for(i in 1:5){
aa <- append(aa, c(df[i,1]:df[i,2]))
i+1
}

hist(aa)
dens <- density(aa)
plot(dens)
$\endgroup$
3
  • $\begingroup$ Thanks @b.nota I am giving this a shot right now and will see how it turns out. This may be exactly what I am looking for. $\endgroup$ – bluescholar1212 Jul 27 '17 at 16:33
  • $\begingroup$ I'm currently working on making the density plot and having some issues. I have a new column with "qstart:end" for each row in my data frame, but when I try to make a density plot with this column R throws an error. How can I get around this? $\endgroup$ – bluescholar1212 Jul 27 '17 at 18:21
  • $\begingroup$ Great man! There are many ways in R how to get the data and make a plot from it. I'll give you a small example in the edit of my answer. Good luck! $\endgroup$ – benn Jul 28 '17 at 19:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.