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

<p>Something like this:</p>

<pre><code>blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'
</code></pre>

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

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

<p>Let me show an example with a small data frame:</p>

<pre><code>qstart &lt;- c(200, 300, 250, 400, 500)
qend &lt;- c(300, 450, 400, 600, 650)

df &lt;- as.data.frame(cbind(qstart,qend))

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

hist(aa)
dens &lt;- density(aa)
plot(dens)
</code></pre>

CC BY-SA 3.0
2017-07-27T16:33:11.803Z
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.
bluescholar1212
CC BY-SA 3.0
2017-07-27T18:21:55.827Z
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?
bluescholar1212
CC BY-SA 3.0
2017-07-28T19:19:37.610Z
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!
benn
CC BY-SA 3.0
2017-07-27T08:55:09.647Z
2017-07-28T19:18:37.343Z
2017-07-28T19:18:37.343Z
benn
<p>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.</p>

<p>For example, from my BLAST search:</p>

<p><a href="https://i.sstatic.net/zVWjr.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zVWjr.png" alt="enter image description here"></a></p>

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

<p>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:</p>

<p><a href="https://i.sstatic.net/HhPYg.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/HhPYg.png" alt="enter image description here"></a></p>

<p>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.</p>

<p>Thanks.</p>

CC BY-SA 3.0
2017-07-27T03:53:28.867Z
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](https://en.wikipedia.org/wiki/List_of_alignment_visualization_software) help
l0o0
CC BY-SA 3.0
2017-07-27T09:06:24.047Z
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?
terdon
CC BY-SA 3.0
2017-07-27T16:32:25.373Z
@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](https://bcbio.wordpress.com/2009/02/07/automated-protein-conservation-display-from-blast-alignments/) provides a description of the method used to make the conservation plot in my post in detail.
bluescholar1212
CC BY-SA 3.0
2017-07-27T16:35:58.313Z
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.
terdon
CC BY-SA 3.0
2017-07-27T16:41:10.427Z
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?
bluescholar1212
<r><blast><visualization>
creating graph of distribution of blast hits on the query sequence
# creating graph of distribution of blast hits on the query sequence

(Asked by: username_0 on 2017-07-26)

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](https://i.sstatic.net/zVWjr.png)](https://i.sstatic.net/zVWjr.png)

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](https://i.sstatic.net/HhPYg.png)](https://i.sstatic.net/HhPYg.png)

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.

**Tags:** r, blast, visualization
> **Comments:**
> **username_1**: 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](https://en.wikipedia.org/wiki/List_of_alignment_visualization_software) help  
> **username_2**: 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?  
> **username_0**: @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](https://bcbio.wordpress.com/2009/02/07/automated-protein-conservation-display-from-blast-alignments/) provides a description of the method used to make the conservation plot in my post in detail.  
> **username_2**: 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.  
> **username_0**: 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?  

---

## 1 answers

### Answer by username_3 on 2017-07-27. Score: 4 | Accepted answer

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)

> **Comments:**
> **username_0**: 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.  
> **username_0**: 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?  
> **username_3**: 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!  

---

