<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/how-to-get-the-homology-of-a-antibody-using-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://r-posts.com/how-to-get-the-homology-of-a-antibody-using-r/"> R-posts.com</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<div>
<div>
<h3>What is homology?</h3>
<p><span>P</span>roteins are conserved bio molecules present in all organisms. Some proteins are conserved in many similar species, making them homologues of each other, meaning that the sequence is not the same but there is a degree of similarity.</p>
<p>This is measured as homology in percentage.</p>
<figure>
<div>
<div></div>
<div><img src="https://i1.wp.com/cdn-images-1.medium.com/max/800/1*F_Ph5K6_6dG_ZX2KO1_qiw.png?w=660&ssl=1" data-recalc-dims="1"></div>
</div><figcaption>Histone H1 homology in mammals (<a rel="nofollow" href="https://en.wikipedia.org/wiki/Sequence_homology" target="_blank">https://en.wikipedia.org/wiki/Sequence_homology</a>)</figcaption></figure>
<h4>What is R</h4>
<p>R is a high level programming language with a focus on statistical analysis and machine learning.<span> </span><a rel="nofollow" href="https://www.r-project.org/about.html" target="_blank">R Project</a></p>
<p>R is widely used in bioinformatics and has a huge resources library of packages, most of them developed by the top Universities to help with hard to understand data, like genomic and proteomic data.</p>
<h4>Why this matters?</h4>
<blockquote><p>This matters because antibodies will bind to specific proteins, but many proteins have similar structure and sequence. So when performing an assay you want to know what specific protein the antibody binds to.</p></blockquote>
<blockquote><p>This is also relevant when choosing an antibody that works in many species, specially when doing research in less used species, like veterinary studies.</p></blockquote>
<h3>Getting Started</h3>
<p>You need the sequence of the epitope used to produce the antibody, I need to know where this product will bind in the target protein. If the sequence is very conserved in mammalians for example, this will be homologue in other species; eg. actins, membrane proteins, etc.</p>
<h4>Obtaining the data</h4>
<p>Data for sequence can be provided by the manufacturer, in some cases it will be in the terms of region of the protein, eg. 162–199aa in C-terminal of RAB22A human.</p>
<p>In these cases I will need to retrieve the protein sequence from<span> </span><a rel="nofollow" href="https://www.uniprot.org/" target="_blank">UNIPROT</a><span> </span>and then cut it according to the information I have on the immunogen, this can be done using a program sequentially, like R<span> </span><a rel="nofollow" href="http://www.gagolewski.com/software/stringi/" target="_blank">stringi</a><span> </span>package or using an excel sheet with<span> </span><em>MID</em><span> </span>function.</p>
<h4>Prepare environment in R (libraries)</h4>
<p>I will use <a rel="nofollow" href="https://www.bioconductor.org/" target="_blank">Bioconductor</a> packages:</p>
<hr>
<pre name="”d26d”">#Install packages: rBLAST, Biostrings, annotate 
install.packages(&quot;devtools&quot;)
library(devtools)
install_github(&quot;mhahsler/rBLAST&quot;)
source(&quot;https://bioconductor.org/biocLite.R&quot;)
biocLite(&quot;Biostrings&quot;)
biocLite(&quot;annotate&quot;)</pre>
<p>These libraries will help us retrieve the homology data, and then I will use simple commands to clean and present the data in a nice way.</p>
<h4>Cleaning data</h4>
<p>Cleaning the data is usually 80% of the work of any data scientist.</p>
<p>I want the homology in simple terms of<span> </span><strong>aa</strong>, or in bioinformatic terms FASTA; eg.</p>
<pre name="”f32f”">ATTAAGCTTATTGC</pre>
<p>You can have terminal places, that need to be removed, usually in terms of ‘-C’, I can remove this using a simple<span> </span><em>gsub</em><span> </span>function.</p>
<pre name="”a146″">dataframe$sequenceImunogen &lt;- gsub(&quot;\\-C&quot;, &quot;&quot;, dataframe$sequenceImunogen)</pre>
<p>For this case the sequence is in a column called<span> </span><em>sequenceImunogen</em><span> </span>that is why I have to make a loop so the protein blast goes through all the sequences in the data frame.</p>
<h4>The function</h4>
<p>I will use the function<span> </span><strong>blastSequences</strong><span> </span>to blast the protein and obtain the homology of this specific sequence to other species. The parameters of the function can be found in the<span> </span><a rel="nofollow" href="https://github.com/mhahsler/rBLAST" target="_blank">literature</a>. I used<span> </span><em>program = “blastp”</em><span> </span>so I perform a<span> </span><strong>P</strong>rotein Blast, I use a hitListSize of 50 and a timeout of 200 seconds, this is from previous experience, as I will do this for multiple sequences I prefer the data to be<span> </span><em>as = data.frame</em>.</p>
<pre name="”420a”">x &lt;- blastSequences(x = Sequence, timeout = 200, as = &quot;data.frame&quot;, program = &quot;blastp&quot;, hitListSize = 50)</pre>
<h4>Making a loop in R</h4>
<p>When creating a loop in R I always put a<span> </span><em>Start</em><span> </span>variable with the<span> </span><em>Sys.time()</em>constant, so I know how long it takes to run this loop.</p>
<pre name="”8c46″">Start &lt;- Sys.time()</pre>
<p>I have use this in a couple of hundred sequences and the result toke a few days, so keep that in mind when you run.</p>
<p>I first create a function that formats my results in terms of percentage of homology and give it the term percentage.</p>
<pre name="”e90e”">#Create a percent function
percent &lt;- function(x, digits = 1, format = &quot;f&quot;, ...) {
paste0(formatC(100 * x, format = format, digits = digits, ...), &quot;%&quot;)
           }</pre>
<h4>Cleaning the result data</h4>
<p>Once the blastP results come, the 1st step is to clean unwanted information so it is not confusing, I used the base functions<span> </span><em>subset</em>,<span> </span><em>grepl</em>,<span> </span><em>gsub</em><span> </span>to get only information with results verified, and remove annoying characters. The objective is to have a nice table result with the best homologues.</p>
<pre name="”ffe1″">#Remove what I do not want
x &lt;- subset(x, grepl(&quot;_A&quot;, x$Hit_accession) == FALSE)
x &lt;- subset(x, grepl(&quot;XP_&quot;, x$Hit_accession) == FALSE)
x &lt;- subset(x, grepl(&quot;unnamed&quot;, x$Hit_def) == FALSE)

#Removing characters that do not 
x$Species &lt;- sub(&quot;.*\\[(.*)\\].*&quot;, &quot;\\1&quot;, x$Hit_def, perl=TRUE)
x$Species &lt;- gsub(&quot;/].*&quot;,&quot;&quot;, x$Species)
x$Hit_def &lt;- gsub(&quot;[[:punct:]].*&quot;, &quot;&quot;, x$Hit_def)</pre>
<h4>Retrieving the results above 80% identical</h4>
<p>After I get all the results and clean a bit of the results, I want to subset the homologues to above 80%. To do this I count the total characters of the sequence using<span> </span><em>nchar</em><span> </span>function and divide the total number of hit results by this number of character, ie. for a 100 aa sequence if I hit 80 I get a 80% homology.</p>
<pre name="”d5e4″">x$Hsp_identity &lt;- as.numeric(x$Hsp_identity)
n &lt;- nchar(dataframe$sequenceImunogen[i])
x$percentage &lt;- as.numeric(x$Hsp_identity/n)
x&lt;- x[order(-x$percentage), ]
x&lt;- x[x$percentage > 0.8, ]
x$percentage &lt;- percent(x$percentage)</pre>
<p>I then clean the result to show as percentage.</p>
<h4>Format nice output</h4>
<p>I want this to be able to be presented in HTML, so I use HTML tags to create a visually pleasing table.</p>
<p>Because this is to be online I want to add hyperlinks to the results so my customers can see the result and analyse the data.</p>
<pre name="”1dfc”">x$link &lt;- paste0(&quot;https://www.ncbi.nlm.nih.gov/protein/&quot;,x$Hit_accession)
    x$final &lt;- paste(x$Hit_def, x$Species, x$percentage, sep = &quot; - &quot;)
    x &lt;- subset(x,!duplicated(x$final))
    x$TOTAL &lt;- paste0(&quot;&lt;a href='&quot;,x$link,&quot;'target='_blank'>&quot;,x$final,&quot;&lt;/a>&quot;)</pre>
<p>I Remove also any duplicated information from our table, and paste all the result into a new column, with all results separated by<span> </span><strong><br>
</strong>(line break in HTML), to do this I use the function<span> </span><em>paste</em>, with the option<span> </span><em>collapse = “<br>”</em></p>
<p>This will create hopefully a table like following:</p>
<figure>
<div>
<div></div>
<div><img src="https://i2.wp.com/cdn-images-1.medium.com/max/800/1*7Zqmi3qZrkdKjZfdUpU2zg.png?w=660&ssl=1" data-recalc-dims="1"></div>
</div>
</figure>
</div>
</div>
</section>
<section>
<div>
<hr></div>
<div>
<div>
<h3>All the code</h3>
<p>All the code can be seen bellow: Use<span> </span><em>TryCatch</em><span> </span>for error management, sometimes the functions times out and you need to rerun, you get a prompt asking to skip or run again, a simple Y is enough.</p>
<pre name="”48fa”">Start &lt;- Sys.time()
    for (i in 1:nrow(dataframe)) {
          tryCatch({
          percent &lt;- function(x, digits = 1, format = &quot;f&quot;, ...) { paste0(formatC(100 * x, format = format, digits = digits, ...), &quot;%&quot;)
                    }
            x &lt;- blastSequences(x = dataframe$sequenceImunogen[i], timeout = 200, as = &quot;data.frame&quot;, program = &quot;blastp&quot;, hitListSize = 50)

#Remove what I do not want
      x &lt;- subset(x, grepl(&quot;_A&quot;, x$Hit_accession) == FALSE)
      x &lt;- subset(x, grepl(&quot;XP_&quot;, x$Hit_accession) == FALSE)
      x &lt;- subset(x, grepl(&quot;unnamed&quot;, x$Hit_def) == FALSE)

#Removing characters that do not 
      x$Species &lt;- sub(&quot;.\\[(.)\\].*&quot;, &quot;\\1&quot;, x$Hit_def, perl=TRUE)
      x$Species &lt;- gsub(&quot;/].*&quot;,&quot;&quot;, x$Species)
      x$Hit_def &lt;- gsub(&quot;[[:punct:]].*&quot;, &quot;&quot;, x$Hit_def)
      x$Hsp_identity &lt;- as.numeric(x$Hsp_identity)
      n &lt;- nchar(masterir2$sequenceImunogen[i])
      x$percentage &lt;- as.numeric(x$Hsp_identity/n)
      x&lt;- x[order(-x$percentage), ]
      x&lt;- x[x$percentage > 0.8, ]
      x$percentage &lt;- percent(x$percentage)

      x$link &lt;- paste0(&quot;https://www.ncbi.nlm.nih.gov/protein/&quot;,x$Hit_accession)
      x$final &lt;- paste(x$Hit_def, x$Species, x$percentage, sep = &quot; - &quot;)
      x &lt;- subset(x,!duplicated(x$final))
      x$TOTAL &lt;- paste0(&quot;&lt;a href='&quot;,x$link,&quot;'target='_blank'>&quot;,x$final,&quot;&lt;/a>&quot;)

  dataframe$Homology[i] &lt;- paste(unlist(x$TOTAL), collapse =&quot;&lt;/br>&quot;)                
}, error=function(e{cat(&quot;ERROR :&quot;,conditionMessage(e), &quot;\n&quot;)})
    }
end &lt;- Sys.time()
print(end - Start)</pre>
<h3>Conclusion</h3>
<p>I can use R and the Bioinformatics packages in Bioconductor to retrieve the homology of known epitopes, this help us provide more valuable information to researchers.</p>
<blockquote><p>Homology is important in the niche world of antibodies, there is a limited amount of money for validation and a limited amount of species that the antibodies can be tested on. With this simple program a homology can be check in a click of a button, instead of using</p></blockquote>
<h3>What next?</h3>
<p>This can be incorporated into a small web application that can provide the homology for a specific selected product with a click off a button, saving time to search and blast on the NCBI website.</p>
</div>
</div>
</section>
<section>
<div>
<hr></div>
<div>
<div>
<p>I am very new to the world of R and coding, so any comment is appreciated.</p>
</div>
</div>
</section>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://r-posts.com/how-to-get-the-homology-of-a-antibody-using-r/"> R-posts.com</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/pWsbtT6oXvE" height="1"></a></div>