<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/proteomics-data-analysis-1-3-data-acquisition-and-cleaning/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/proteomics-data-analysis-1-3-data-acquisition-and-cleaning/"> R Programming – DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>The analysis of DNA and RNA, the blueprint of life and its carbon copy, has become a staple in the burgeoning field of molecular biology. An emerging and exciting area of study that adds another dimension to our understanding of cellular biology is that of proteomics. The use of mass spectrometry has enabled the identification and quantification of thousands of proteins in a single sample run. </p>
<p>In this tutorial series, I will break down the steps of processing a high-throughput proteomics data set derived from mass spectrometry analysis as follows: </p>
<ul>
<li>Data acquisition and cleaning</li>
<li>Data filtering and missing value imputation</li>
<li>Data visualization and interpretation</li>
</ul>
<h2>Source of Proteomics Data</h2>
<p>To obtain a sample data set, I combed through a proteomics data repository called <a rel="nofollow" href="https://www.ebi.ac.uk/pride/archive/simpleSearch" target="_blank">PRIDE</a> and found an interesting study on drug resistance in breast cancer cell lines. I downloaded the <a rel="nofollow" href="https://www.ebi.ac.uk/pride/archive/projects/PXD002057/files" target="_blank">raw files</a>, which are the output from mass spectrometry analysis, and processed them using a software called <a rel="nofollow" href="http://www.biochem.mpg.de/5111795/maxquant" target="_blank">MaxQuant</a> to map the spectral data to protein sequences. A total of six raw files, corresponding to three replicates of two cell lines (one resistant, one non-resistant), was used. There are numerous other ways to process mass spectrometry data (e.g. Mascot, SEQUEST, ProteinProspector) so the final data table of protein abundance measurements will vary base on the approach. The starting point for this tutorial is the MaxQuant <em>ProteinGroups</em> output file, which can be downloaded <a rel="nofollow" href="https://ucsf.box.com/s/1wt4hjdt3yq9g8fcs8oam1btg66dzt54" target="_blank">here</a>.</p>
<h2>Data Acquisition</h2>
<p>The first step is to read the tab-separated data file into R.</p>
<pre><code>raw = read.delim(&quot;proteinGroups.txt&quot;, stringsAsFactors = FALSE, colClasses = &quot;character&quot;)
</code></pre>
<p>Our raw data is an enormous 1787-by-79 data frame. Proteins are arranged in rows and their information in columns. The primary columns of interest are those containing intensity measurements, which reflects protein abundance.</p>
<pre><code>grep(&quot;^LFQ.intensity&quot;, names(raw), value = TRUE)</code>
<em>## [1] &quot;LFQ.intensity.Parental_bR1&quot; &quot;LFQ.intensity.Parental_bR2&quot; ## [3] &quot;LFQ.intensity.Parental_bR3&quot; &quot;LFQ.intensity.Resistant_bR1&quot;
## [5] &quot;LFQ.intensity.Resistant_bR2&quot; &quot;LFQ.intensity.Resistant_bR3&quot;
</em></pre>
<p>Again, we have a total of six samples. The <em>Parental</em> represents intensity data from the breast cancer cell line SKBR3 while the <em>Resistant</em> is an drug-resistant cell line derived from culturing the parentals in the presence of an EGFR inhibitor. For more information regarding the study, please see the original <a rel="nofollow" href="http://europepmc.org/abstract/MED/26883193" target="_blank">publication</a>.</p>
<h2>Data Clean Up</h2>
<h3>Remove False Hits</h3>
<p>The next step after data acquisition is to clean and organize our data. The first order of business is to remove the false hits, including the contaminants, reverse proteins, and proteins identified by site. These are annotated with a “+” under the columns <em>Potential.contaminant</em>, <em>Reverse</em>, and <em>Only.identified.by.site</em>. We filter the data frame by keeping rows without a “+” annotation in any of the three columns.</p>
<pre><code>library(dplyr)
df = dplyr::filter(raw, Potential.contaminant != &quot;+&quot;) %>%
  dplyr::filter(Reverse != &quot;+&quot;) %>%
  dplyr::filter(Only.identified.by.site != &quot;+&quot;)
</code></pre>
<p>Often, there is a column indicating the confidence of the protein identification. In our case, that is <em>Q.value</em>, which represents the probability that the protein is a false hit. A typical cutoff is set at 0.01. Fortunately, MaxQuant takes care of this operation and all Q values are below the threshold.</p>
<pre><code>summary(as.numeric(df$Q.value))</code>
<em>## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0000000 0.0000000 0.0000000 0.0007193 0.0012759 0.0091429
</em></pre>
<h3>Define Protein and Gene IDs</h3>
<p>A quick look at <em>Protein.IDs</em> and <em>Fasta.headers</em> columns tells us that the protein IDs, protein names, and gene IDs are not clearly displayed.</p>
<pre><code>head(df$Protein.IDs)</code>
<em>## [1] &quot;sp|A0AV96|RBM47_HUMAN&quot; &quot;sp|A0FGR8|ESYT2_HUMAN&quot; &quot;sp|A1L0T0|ILVBL_HUMAN&quot;
## [4] &quot;sp|A4D1S0|KLRG2_HUMAN&quot; &quot;sp|A5YKK6|CNOT1_HUMAN&quot; &quot;sp|A6NDG6|PGP_HUMAN&quot;
</em></pre>
<pre><code>head(df$Fasta.headers)</code>
<em>## [1] &quot;>sp|A0AV96|RBM47_HUMAN RNA-binding protein 47 OS=Homo sapiens GN=RBM47 PE=1 SV=2&quot; ## [2] &quot;>sp|A0FGR8|ESYT2_HUMAN Extended synaptotagmin-2 OS=Homo sapiens GN=ESYT2 PE=1 SV=1&quot; ## [3] &quot;>sp|A1L0T0|ILVBL_HUMAN Acetolactate synthase-like protein OS=Homo sapiens GN=ILVBL PE=1 SV=2&quot; ## [4] &quot;>sp|A4D1S0|KLRG2_HUMAN Killer cell lectin-like receptor subfamily G member 2 OS=Homo sapiens GN=KLRG2 PE=1 SV=3&quot;
## [5] &quot;>sp|A5YKK6|CNOT1_HUMAN CCR4-NOT transcription complex subunit 1 OS=Homo sapiens GN=CNOT1 PE=1 SV=2&quot; ## [6] &quot;>sp|A6NDG6|PGP_HUMAN Glycerol-3-phosphate phosphatase OS=Homo sapiens GN=PGP PE=1 SV=1&quot;
</em></pre>
<p>We will use regular expressions to extract the protein names into a column named <em>Protein.name</em>, UniProt protein IDs into <em>Protein</em>, and gene IDs into <em>Gene</em>. Note that some rows may have multiple IDs separated by semicolons. In those instances, we will isolate the first one available.</p>
<pre><code># Isolate the first ID
df$Protein.IDs = sub(&quot;;.*&quot;, &quot;&quot;, df$Protein.IDs)
df$Fasta.headers = sub(&quot;;.*&quot;, &quot;&quot;, df$Fasta.headers)
# Extract Protein name
regex = regexpr(&quot;(?&lt;=_HUMAN.).*(?=.OS)&quot;, df$Fasta.headers, perl = TRUE)
df$Protein.name = regmatches(df$Fasta.headers, regex)
# Extract UniProtID
regex = regexpr(&quot;(?&lt;=\\|).*(?=\\|)&quot;, df$Protein.IDs, perl = TRUE)
df$Protein = regmatches(df$Protein.IDs, regex)
# Extract Gene ID
regex = regexpr(&quot;((?&lt;=\\|[[:alnum:]]{6}\\|).*(?=_HUMAN)|(?&lt;=\\|[[:alnum:]]{10}\\|).*(?=_HUMAN))&quot;,
                df$Protein.IDs, perl = TRUE)
df$Gene = regmatches(df$Protein.IDs, regex)
</code></pre>
<h3>Transform Intensity Columns</h3>
<p>Due to our function call for reading the data table, all columns are cast as the <strong>character</strong> data type. It will benefit our analysis to convert the intensity columns to the proper <strong>numeric</strong> data type.</p>
<pre><code>intensity.names = grep(&quot;^LFQ.intensity&quot;, names(raw), value = TRUE)
df[intensity.names] = sapply(df[intensity.names], as.numeric)
</code></pre>
<p>Now, let us examine the distribution of protein intensities in a sample. Below is a histogram of the <em>Parental_bR1</em> protein intensities.</p>
<p><a rel="nofollow" href="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/08/unnamed-chunk-9-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/08/unnamed-chunk-9-1-490x490.png?w=450&ssl=1" alt="" class="wp-image-16240" data-recalc-dims="1"></a></p>
<p>The distribution is clearly skewed to the right with a few highly abundant proteins. To normalize the distribution, it is common practice to log<sub>2</sub> transform the intensity data.</p>
<pre><code>LOG.names = sub(&quot;^LFQ.intensity&quot;, &quot;LOG2&quot;, intensity.names)   # rename intensity columns
df[LOG.names] = log2(df[intensity.names])
</code></pre>
<p>Let's check out the transformed distribution on <em>Parental_bR1</em> (much better!):</p>
<p><a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/08/unnamed-chunk-11-1.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/08/unnamed-chunk-11-1-490x490.png?w=450&ssl=1" alt="" class="wp-image-16241" data-recalc-dims="1"></a></p>
<h2>Summary</h2>
<p>This is the first of three tutorials on proteomics data analysis. I have outlined the steps to read and clean a typical mass spectrometry-based proteomics data set. </p>
<p>In the next tutorial, we will examine the data in greater detail. In doing so, we will find that only a handful of proteins are quantified across all samples. In other words, proteins are often picked up in one sample but not in the others. This is known as the missing value problem. Stick around to explore techniques for filtering proteins based on the number of valid values and filling in the missing values using data imputation.</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p></p></strong>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/accessing-web-data-json-in-r-using-httr/" title="Permanent Link toAccessing Web Data (JSON) in R using httr" target="_blank">Accessing Web Data (JSON) in R using httr</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/zomato-web-scraping-with-beautifulsoup-in-python/" title="Permanent Link toZomato Web Scraping with BeautifulSoup in Python" target="_blank">Zomato Web Scraping with BeautifulSoup in Python</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/processing-huge-dataset-with-python/" title="Permanent Link toProcessing Huge Dataset with Python" target="_blank">Processing Huge Dataset with Python</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/finding-the-popular-indian-blogging-platform-by-web-scraping-in-python/" title="Permanent Link toFinding the Popular Indian Blogging Platform by Web Scraping in Python" target="_blank">Finding the Popular Indian Blogging Platform by Web Scraping in Python</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/osm-nominatim-with-r-getting-locations-geo-coordinates-by-its-address/" title="Permanent Link toOSM Nominatim with R: getting Location’s Geo-coordinates by its Address" target="_blank">OSM Nominatim with R: getting Location’s Geo-coordinates by its Address</a></li>
</ol>
	<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://datascienceplus.com/proteomics-data-analysis-1-3-data-acquisition-and-cleaning/"> R Programming – DataScience+</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/qGonZBdHFC0" height="1"></a></div>