# creating graph of distribution of blast hits on the query sequence

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:

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:

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.

• 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
– l0o0
Jul 27, 2017 at 3:53
• 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? Jul 27, 2017 at 9:06
• @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. Jul 27, 2017 at 16:32
• 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. Jul 27, 2017 at 16:35
• 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? Jul 27, 2017 at 16:41

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)

• 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. Jul 27, 2017 at 16:33
• 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? Jul 27, 2017 at 18:21
• 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
Jul 28, 2017 at 19:19