Received on XXXXX; revised on XXXXX; accepted on XXXXX

Associate Editor: XXXXXXX

# Mass spectrometry based protein identification with accurate statistical
significance assignment

Gelio Alves and Yi-Kuo Yu111to whom correspondence should be addressed
National Center for Biotechnology Information,
National Library of Medicine, National Institutes of Health
8600 Rockville Pike, Bethesda, MD 20894, USA

(2014; 2014)

###### Abstract

##  1 Motivation:

Assigning statistical significance accurately has become increasingly
important as meta data of many types, often assembled in hierarchies, are
constructed and combined for further biological analyses. Statistical
inaccuracy of meta data at any level may propagate to downstream analyses,
undermining the validity of scientific conclusions thus drawn. From the
perspective of mass spectrometry based proteomics, even though accurate
statistics for peptide identification can now be achieved, accurate protein
level statistics remain challenging.

##  2 Results:

We have constructed a protein ID method that combines peptide evidences of a
candidate protein based on a rigorous formula derived earlier; in this formula
the database $P$-value of every peptide is weighted, prior to the final
combination, according to the number of proteins it maps to. We have also
shown that this protein ID method provides accurate protein level $E$-value,
eliminating the need of using empirical post-processing methods for type-I
error control. Using a known protein mixture, we find that this protein ID
method, when combined with the Sorić formula, yields accurate values for the
proportion of false discoveries. In terms of retrieval efficacy, the results
from our method are comparable with other methods tested.

##  3 Availability:

The source code, implemented in C++ on a linux system, is available for
download at ftp://ftp.ncbi.nlm.nih.gov
/pub/qmbp/qmbp_ms/RAId/RAId_Linux_64Bit

##  4 Contact:

yyu<EMAIL_ADDRESS>
##  5 Introduction

Peptide identifications (ID) via mass spectrometry (MS) have become the
central component in modern proteomics; this component, combined with
additional analyses, routinely yields pragmatic meta data, including protein
ID, protein quantification, protein structure and protein associations (Zhang
et al., 2013). These meta data, especially the associated statistical
significance assignments, need to be as accurate as possible because they
often form the building blocks for investigations at the systems biology level
and influence the scientific conclusions drawn henceforth. In this paper, we
focus on protein ID, in particular on improving the accuracy of statistical
significance assigned to proteins identified.

The need for robust developments towards accurate statistical significance
assignments has been advocated (Noble and MacCoss, 2012; Huang et al., 2012)
despite the existence of many protein ID methods (McHugh and Arthur, 2008;
Serang and Noble, 2012; Li and Radivojac, 2012). It has also been suggested
(Spirin et al., 2011) that the primary cause of unreliable significance
assignment for protein ID can be attributed to inaccurate significance
assignment for peptide ID. Frequently used error-control/significance-
assigning methods for peptide ID largely fall into two groups: proportion of
false discovery (PFD), which is often incorrectly termed as false discovery
rate (Benjamini and Hochberg, 1995), and spectrum-specific $P$-value/$E$-value
(Fenyo and Beavis, 2003; Alves et al., 2007; Park et al., 2008). Methods
belonging to the first group, controlling type-I error globally only, do not
discriminate among identified peptides (Elias and Gygi, 2007). Methods
belonging to the second group, capable of assigning per-spectrum per-peptide
significance, can properly prioritize identified peptides when reported
$P$-/$E$-values are accurate; but the needed statistical accuracy is often
unattainable due to improper heuristics or unjustifiable distribution
assumptions (Alves et al., 2007; Segal, 2008; Spirin et al., 2011).

Given a tandem MS (MS/MS) spectrum and a quality score cutoff ${\sf S}_{c}$,
the $E$-value $E({\sf S}_{c})$ should reflect the expected number of random
peptides with scores the same as or better than ${\sf S}_{c}$. (Similarly, the
$P$-value $P({\sf S}_{c})$ reflects the probability of finding a random
peptide with quality score ${\sf S}\geq{\sf S}_{c}$.) In general, the
$E$-value is obtained by multiplying the $P$-value by the total number of
qualified peptides (whose masses fall in the range
$[m_{p}-\delta,m_{p}+\delta]$ with $m_{p}$ being the precursor ion’s mass and
$\delta$ the specified tolerance) in the database searched. Thus, besides
providing the user with the numbers of false positives to anticipate, accurate
$E$-value assignments enable ranking of candidate peptides across different
spectra and experiments. In database searches in proteomics, the goal of
accurate statistics can be approached in at least two ways. First, one may
devise a scoring function whose resulting score distribution can be
analytically characterized and thus used to infer the statistical significance
(Alves et al., 2007); if this is done correctly, the theoretical score
distribution should fit well the bulk part of the normalized score histogram
obtained from scoring all qualified peptides in the database of interest.
Second, one may infer the spectrum-specific $P$-value via the normalized score
histogram obtained from scoring all possible peptides (APP) (Alves and Yu,
2008); in this case, the database dependence appears only in the $E$-value,
which is the $P$-value multiplied by the number of qualified peptides
associated with the specified precursor ion mass and mass error tolerance.
Either way yields database-specific $E$-values. Once a peptide $E$-value is
obtained, one may transform it into the peptide database $P$-value (DPV) (Yu
et al., 2006; Alves et al., 2008b), representing the likelihood of obtaining,
in the database chosen, at least one peptide scoring equal to or better than
the prescribed threshold. When combining $P$-values of peptides associated
with a candidate protein, we use the peptides’ DPVs.

Specifically, our proposed protein ID method combines peptide evidences of a
candidate protein using a rigorous formula derived earlier (Alves and Yu,
2011); in this formula the DPV of every peptide is weighted, prior to the
final combination, according to the number of proteins it maps to. Among the
existing protein ID methods, the approach taken by Spirin et al. (2011) is
closest to ours; both methods combine peptides’ spectrum-specific $P$-values.
There are, however, major differences between our method and that of Spirin et
al. (2011). First, in our method, each candidate peptide of a query spectrum
receives a DPV, allowing multiple matching peptides per spectrum. This is to
take into account the possibility of peptide co-elution (Alves et al., 2008a).
For the method of Spirin et al. (2011), only the best peptide match per
spectrum is considered and the peptide DPV thus represents the probability of
having the best match score no worse than the prescribed threshold when
searching a database. Since each random protein database only contributes one
best match score, searching many random protein databases is required for the
$P$-value assignment. Second, the candidate peptides’ $P$-values are combined
differently. Our method, down-weighting contributions of peptides mappable to
multiple proteins, combines peptide DPVs directly using a rigorous formula
(Alves and Yu, 2011); the method of Spirin et al. (2011) first transforms, for
every candidate protein, the $P$-values of its associated peptides into
$Z$-scores, combines them using Stouffer’s formula (Whitlock, 2005), and then
transforms the combined $Z$-score back to a final $P$-value with multiple
hypotheses testing correction. Third, the cutoff conditions for peptides’
$P$-values are different. Our method approximates DPVs (Yu et al., 2006; Alves
et al., 2008b) by $E$-values, valid for small $E$-values, and retains all
peptides whose $E$-values are less than one. That is, we have a global cutoff
condition. For the method of Spirin et al. (2011), the peptide cutoff
$P$-value varies by candidate protein: given a candidate protein, its
corresponding peptides’ $Z$-scores are first sorted in descending order; the
$k$th $Z$-score is chosen as the cutoff provided that the maximum combined
$Z$-score is reached while combining the top $k$ $Z$-scores using the
Stouffer’s formula.

There exist many other protein ID methods, for example, ProFound (Zhang and
Chait, 2000), ProteinProphet (Nesvizhskii et al., 2003), DBParser (Yang et
al., 2004), EBP (Price et al., 2007), PANORAMICS (Feng et al., 2007), PROVALT
(McHugh and Arthur, 2008), X!Tandem (Fenyo et al., 2010), Scaffold (Searle,
2010) and npCI (Serang et al., 2013), to name just a few. We refer the readers
to recent review papers (Huang et al., 2012; Serang and Noble, 2012) for
details and more comprehensive listings of these methods. Although some of
them do start with spectrum-specific peptide $P$-values, they often assume
certain parametric forms for the peptide score distributions when searching a
random database; other methods, however, only process outputs of specific
peptide identification tools, limiting their uses to certain platforms. By
discarding all but the best few peptide scores per spectrum per database
search, the method of Spirin et al. (2011) does not rely on the accuracy of
the full peptide score distribution from searching a random database and in
principle can accept input from various peptide identification tools. Our
method is free from the aforementioned problems for different reasons. Founded
on a derived analytical formula, our method can be applied in general and will
yield accurate protein $P$-values if the input peptide DPVs (or $E$-values)
are accurate. When using peptide $E$-values reported by RAId_DbS, even though
the parameters of the score distribution are determined via maximum-
likelihood, the functional form of the score distribution is analytically
derived (Alves et al., 2007) rather than assumed. When the statistical
significances are obtained from RAId_aPS (Alves et al., 2010), for every
scoring function implemented, the $P$-values are inferred by scoring APP
instead of assuming that the score histogram follows a specific form; the
peptide $E$-values are then obtained via multiplying the $P$-values by the
respective numbers of qualified peptides.

The paper is organized as follows. The mathematical underpinnings of our
formalism will be described in the methods section. In the results section,
comparisons of our method with other approaches will be made; the accuracy of
the reported protein $P$-value will be illustrated. Some technical but
important issues will be addressed in the discussion section. To keep the
paper focused, we relegate to supplementary information figures and tables
that complement or corroborate the information contained in the main text.

##  6 Methods

###  6.1 Statistical Protocols

Weighting the contribution of each peptide in protein identification is
important. It helps mitigate the issue of peptide degeneracy, where an
identified peptide is a subsequence of multiple database proteins. The optimal
weighting scheme, however, can depend on the protein ID methodology employed.
For the purpose of our study, namely, devising a method that yields accurate
protein $P$-values, we opt for a simple weighting scheme: a peptide’s weight
is inversely proportional to the number of database proteins it maps to.
Within a sample, when multiple spectral searches identify the same peptide but
with different significance levels, only the most significant assignment of
that peptide is retained for further analyses.

The foundation of our method is built upon a rigorous formula (Mathai, 1983;
Alves and Yu, 2011) that enables weighted combination of $P$-values. When the
weights are all identical, this formula reduces to Fisher’s formula (Fisher,
1932; Bahrucha-Reid, 1960); when the weights are all different, this formula
reduces to the formula of Good (Good, 1955). A detailed derivation and
generalization to incorporating nearly identical weights can be found in
(Alves and Yu, 2011), whose notation will be used to briefly summarize the
content of the formula.

Let us assume that a given protein contains $L$ identified peptides with
$P$-values. Let us further group these $L$ peptides, according to the number
of database proteins a peptide maps to, into $m$ groups with $1\leq m\leq L$.
Within each group $k$, the $n_{k}$ peptide $P$-values are weighted equally;
while peptide $P$-values in different groups are weighted differently.

The weighting enters our formalism through the following quantities of
interest

$\displaystyle\tau$ $\displaystyle\equiv$
$\displaystyle\prod_{k=1}^{m}\left[\prod_{j=1}^{n_{k}}p_{k;j}\right]^{w_{k}}\;,$
(1) $\displaystyle Q$ $\displaystyle\equiv$
$\displaystyle\prod_{k=1}^{m}\left[\prod_{j=1}^{n_{k}}x_{k;j}\right]^{w_{k}}\;,$
(2)

where each $p_{k;j}$ represents a reported peptide $P$-value, each $x_{k;j}$
represents a random variable drawn from an uniform, independent distribution
over $(0,1]$ and each $w_{k}$ is a positive weight. The quantity of interest
${\rm Prob}(Q\leq\tau)$, representing the protein $P$-value, was obtained
earlier (Alves and Yu, 2011) and is repeated below for clarity.

Let $F(\tau)\equiv{\rm Prob}(Q\leq\tau)$, one may show that

$\displaystyle F(\tau)$ $\displaystyle=$
$\displaystyle\left[\prod_{l=1}^{m}r_{l}^{n_{l}}\right]\sum_{k=1}^{m}\,\sum_{\mathcal{G}(k)}\;\biggr{\\{}\frac{1}{r_{k}^{g_{k}+1}}H(-r_{k}\ln\tau\,,\,g_{k})\times$
(3) $\displaystyle\times\biggr{(}\prod_{j=1,j\neq
k}^{m}\frac{(n_{j}-1+g_{j})!}{(n_{j}-1)!g_{j}!}\frac{(-1)^{g_{j}}}{(r_{j}-r_{k})^{n_{j}+g_{j}}}\biggr{)}\biggr{\\}},$

where $r_{k}\equiv 1/w_{k}$ is the number of proteins a group-$k$ peptide maps
to, $\sum_{\mathcal{G}(k)}$ enumerates each set of nonnegative integers
$\\{g_{1},g_{2},\ldots,g_{m}\\}$ that satisfies the $k$-dependent constraint
$\sum_{i=1}^{m}g_{i}=n_{k}-1$, and the function $H$ is defined as

$H(x,n)\,\equiv\;e^{-x}\,\sum_{k=0}^{n}\frac{x^{k}}{k!}\;.$ (4)

See the supplementary information for an example application of formula (3).

When searching a database with a prescribed peptide mass error tolerance
$\delta$, one often needs to score different numbers of database peptides for
spectra with different precursor ion masses. That is, the number of tested
hypotheses (database peptides in the mass range $[m_{p}-\delta,m_{p}+\delta]$)
varies by the precursor ion mass $m_{p}$. The effect of varying number of
multiple hypotheses tested can be properly accounted for by using the peptide
DPVs (Yu et al., 2006; Alves et al., 2008b) for $P$-values ($p_{k;j}$) in eq.
(1); given a quality score cutoff ${\sf S}_{c}$, the peptide DPV is defined as

$P_{\rm db}({\sf S}_{c})=1-e^{-E({\sf S}_{c})}\;,$ (5)

where $E({\sf S}_{c})$ represents the expected number of peptides having score
${\sf S}\geq{\sf S}_{c}$, and the DPV $P_{\rm db}({\sf S}_{c})$ represents the
probability of seeing one or more peptides in a given random database with
quality scores ${\sf S}\geq{\sf S}_{c}$. Another advantage of using DPV is
that as a function of the quality score ${\sf S}$, the $E$-value $E({\sf S})$,
determined by the search score histogram per spectrum and the number of
qualified peptides (database-dependent), correctly takes into account both the
spectrum-specificity and the database-specificity of scoring statistics.

Since the $E$-value specifies the expected number of random database peptides
having scores equal to or better than the given cut-off, a peptide with
$E$-value larger than one is more likely to be a false positive than a true
positive. For this reason, when constructing the evidence peptide set for
identification of a protein, we only include peptides with $E$-values less
than one. This implies that only peptide DPVs less than $(e-1)/e$ are
considered, leading to a combination of truncated $P$-values. Unfortunately,
combining truncated $P$-values, even though doable, is far more complicated
than using eq. (3). However, two observations simplify the matter. First, it
is evident from eq. (5) that the DPV approaches the $E$-value when the
$E$-value is small. Second, we note that confidently identified proteins must
contain evidence peptides with high identification confidences (or small
$E$-values). Therefore, for practical uses, we may approximate the DPV by its
corresponding $E$-value. Because only $E$-values less than one are considered,
the approximated DPVs (or simply the $E$-values) now encompass the full range
between zero and one. Consequently, it is unnecessary to combine truncated
$P$-values, and the simple formula (3) becomes applicable. The protein
$E$-value is then obtained via multiplying the protein $P$-value by a
Bonferroni correction factor; in this case, the Bonferroni factor is the
number of protein clusters (described below) each having at least one evidence
peptide with $E$-value less than one.

We denote by a protein cluster a group of entangled proteins that share a
substantial portion of evidence peptides. To avoid exaggerating the number of
identified proteins, several existing methods (Huang et al., 2012) report
those entangled proteins as one. Adopting the same idea, we implemented this
strategy via a transitive approach described below. One first sorts the
identified proteins by the number of identified evidence peptides in
descending order and using the rank of a protein in the sorted list as that
protein’s cluster index. Starting with the first protein as the reference
protein, all other lower-ranking proteins sharing at least $95\%$ of evidence
peptides with the first protein will have their cluster indexes changed to
that of the reference protein. One then moves the reference point (from the
first) to the second protein, all other lower-ranking proteins sharing at
least $95\%$ of evidence peptides with the reference protein will have their
cluster indexes changed to that of the reference protein. The reference point
is then moved to the third protein and the process continues till the
reference point moves through all proteins in the list. The protein with most
significant $P$-value within a cluster (containing one or more proteins) is
called the head of that cluster, the other proteins members of that cluster.
An exception to the aforementioned clustering rule, however, is introduced to
appropriately emphasize a protein’s evidence peptides that are not shared by
other proteins. We call evidence peptides of this kind unique peptides to a
protein. When a protein has a unique evidence peptide with $E$-value less than
$10^{-4}$, our method doesn’t allow this protein to be a member protein of any
cluster.

###  6.2 MS/MS Datasets

Sixty-three spectral datasets were categorized into four data groups. See
supplementary tables (Table S1 to S4) for details. Protein mixtures giving
rise to spectral datasets were reduced with iodoacetamide, resulting in the
addition of the carbamidomethyl group (57.07 Da) to cystine residues. Each
protein mixture was further digested with trypsin. Among these spectral
datasets, there are also dataset-specific parameters such as the target
database, the maximum number of missed cleavage sites allowed, the precursor-
ion mass error tolerance and the product-ion mass error tolerance. The
dataset-specific parameters are given in the figure caption to provide more
information underlying the generation of the figures.

For brevity, we shall denote the MS/MS spectra obtained from a sample by SN
followed by its sample index. For example, SN1 denotes the collection of MS/MS
spectra acquired from mixture sample one. The first data group, SN1 through
SN15, contained MS/MS spectra from replicates of different dilutions of
Sigma49, a protein standard mixture composed of 49 know human proteins. The
second data group, SN16 through SN26, was downloaded from the Pacific
Northwest National Laboratory and contained spectra from eleven whole-cell-
lysate samples of protein mixtures of Escherichia coli K-12. The third data
group, SN27 through SN30, consisted of spectra from four in-house whole-cell-
lysate samples of protein mixtures of Escherichia coli K-12. Downloaded from
PeptideAtlas database, the fourth data group (SN31 through SN63) was composed
of spectra from SDS-PAGE protein fractionation extractions of human lung
cells.

###  6.3 Protein Databases and Random Databases

Because protein mixtures from Escherichia coli K-12 and Homo sapiens were
analyzed using their corresponding MS/MS spectra, protein databases for both
organisms were thus required. From UniProt
http://www.uniprot.org/downloadshttp://www.uniprot.org/downloads, we
downloaded 4,303 non-redundant protein sequences of Escherichia coli K-12. A
non-redundant Homo sapiens protein database, containing $31,236$ protein
sequences, was obtained from the NCBI site
ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/mRNA_Prot/ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/mRNA_Prot/.

When analyzing statistical significance, it is often required to have random
(decoy) databases in addition to the organismal (target) databases. One common
problem when using random databases is that for a given precursor ion mass the
numbers of qualified peptides in the random database and in the organismal
database may significantly differ. This causes an additional uncertainty in
assessing statistical significance (Elias and Gygi, 2007; Wang et al., 2009).
We can avoid this problem by ensuring that the numbers of qualified peptides
per spectrum are identical for both the random and the organismal databases:
for each qualified peptide in the organismal database, we generate a
corresponding random peptide by randomly shuffling its amino acids.

##  7 Results

The results will be described in the following order. First, we illustrate
that our $E$-value assignments are accurate at both the peptide and the
protein levels. We further show that using the formula proposed by Sorić
(1989), our reported PFDs agree well with the target-decoy PFDs. Second, our
protein $E$-value accuracy is compared with that of using the formulas in
(Spirin et al., 2011). By extending the formula of Sorić for the method of
Spirin et al. (2011), we also evaluate the agreement between their reported
PFDs and the target-decoy PFDs. Benchmarking with some of the existing protein
ID methods will be described in the third part.

###  7.1 $E$-value accuracy

The input peptide DPVs for our protein ID method are obtained via eq. (5)
using the $E$-values reported by RAId_DbS. For this reason, the input peptide
DPVs (for protein ID) are synonymous with the reported peptide DPVs (from
RAId_DbS). As mentioned earlier, the statistical accuracy of our protein ID
method relies on the DPVs for the evidence peptides being accurate. We
therefore start by comparing the input peptide DPV with its definition. In
panel A of Figure 1, the abscissa records the peptide DPV, while the ordinate
displays the observed DPV (i.e., fraction of spectra having at least one or
more matching peptides with reported DPVs smaller than the specified
threshold). The agreement between the observed DPV and the reported DPV
indicates that the peptide DPVs used as input for our protein ID method are
accurate.

Figure 1:  Assessment of $E$-value accuracy. In panels A, B, and C, the closer
the displayed curves are to the $y=x$ line the better. In panel D, the closer
the two displayed curves are to each other the better. See section 7.1 for
more details. For panels A, B, and C, spectral dataset SN 26 (Escherichia coli
K-12 whole cell lysate) is used to search the Escherichia coli database with
mass accuracy $\pm$ 0.033 Da. for precursor and product ions. For panel D,
spectral datasets SN 13-15 (Sigma49 protein standard mixture) is used to
search the Homo sapiens database with precursor ion accuracy $\pm$ 0.033 Da.
and product ion accuracy $\pm$ 0.8 Da.

To assess whether approximating peptide DPVs by their corresponding $E$-values
for $E$-values less than one is reasonable or not, we plot in panel B of
Figure 1 the observed peptide DPVs versus $E$-values. As expected, when
$E$-values are close to one, there is certain degree of disagreement; while
for small $E$-values, the agreement is excellent. To assess the accuracy of
the protein $P$-values reported by our eq. (3), we compare them with the
observed protein $P$-values. As described in the method section, the reported
proteins appear in clusters, each represented by a head protein and its
$P$-value. The observed protein $P$-value is defined as the fraction of
identified protein clusters (whose member proteins each containing at least
one evidence peptide with $E$-value less than one) that have reported
$P$-values smaller than a given threshold. As shown in panel C of Figure 1,
good agreement between the reported protein $P$-value and the observed protein
$P$-value is obtained, indicating that our reported protein $P$-values are
accurate. More protein $P$-value accuracy assessment examples can be found in
Supplementary Figures S1, S2, and S3. With an accurate protein $P$-value, one
can also obtain its corresponding protein $E$-value by multiplying it by the
total number of protein clusters. In Supplementary Figure S4, we show that
reported protein $E$-values obtained this way are accurate.

By having accurate protein $E$-values, one can avoid the uncertainty
associated with using a decoy database (Gupta et al., 2011) while estimating
the proportion of false discoveries. In panel D of Figure 1, we plot two PFD
curves: one is computed using the reported protein $E$-value to estimate the
number of false identifications (hence the PFD), while the other is computed
using the observed PFD obtained from known target protein content in the
sample (Sigma49). The excellent agreement between the observed PFD and the
reported PFD indicates that one should be able to trust the PFD estimated from
accurate reported protein $E$-values. More accuracy assessment examples of the
reported PFD can be found in Supplementary Figure S5.

###  7.2 Comparison with an EVD-based method

Since the method of Spirin et al. (2011) is closest to ours, we also
implemented their method and compute equivalent quantities for comparison.
Following the Supplementary Material of (Spirin et al., 2011), we have
implemented 100 random databases each containing 10,000 random amino acid
sequences. However, instead of generating sequences of uneven length, we opt
for uniform length (each sequence is of length 350) and generate these random
sequences using the background amino acid frequencies of Robinson and Robinson
(1991). The EVD parameters are obtained by using only the best score per
database search and by applying standard procedures described in (Spirin et
al., 2011). The effect of database size difference, leading to rescaling of
the $\alpha$ parameter, is done the same way as in (Spirin et al., 2011).

Figure 2:  Statistical Accuracy Comparison. Except that results from two
methods are being displayed, panels A, B, and C display similar information
respectively to panels A, C and D of Fig. 1. See text in section 7.2 for more
details.

A moment of reflection reveals that the best match $P$-value of (Spirin et
al., 2011) is in fact the DPV (Yu et al., 2006; Alves et al., 2008b). We
therefore plot in panel A of Figure 2 the reported peptide DPVs against the
observed peptide DPVs. The result indicates that the peptide DPV reported by
Spirin et al. (2011) is quite accurate, with an uncertainty of a factor of $2$
as reported by Spirin et al. (2011).

To have a fair assessment, the same procedure for clustering proteins is also
applied to the proteins identified using protocols of Spirin et al. (2011).
The observed protein $P$-value is defined similarly. Database proteins that
contain any of the best match peptides, one from each spectrum, form the
effective protein set. The observed protein $P$-value is simply the fraction
of proteins in the effective protein set that have reported protein $P$-values
less than the specified threshold. The reported protein $P$-value for the head
protein of each cluster is obtained by applying the iterative procedure
(involving uses of Stouffer’s formula) described in (Spirin et al., 2011). In
panel B of Figure 2, the reported protein $P$-values are plotted against the
observed protein $P$-values. The agreement between the reported protein
$P$-values and the observed protein $P$-values is not as great as in the
peptide case. The protein $E$-value is then obtained by multiplying the
protein $P$-value by the total number of proteins in the effective protein
set.

To construct a PFD curve, it is necessary to estimate the number of false
identifications at a given significance threshold. The number of false
identifications can be estimated either by using the reported protein
$E$-values or the number of identifications within the decoy databases. The
latter is currently widely used mainly because accurate protein $E$-values (or
$P$-values) are generally hard to attain. To investigate the agreement between
the PFD curves obtained using decoy databases and using reasonably accurate
protein $P$-values, we use spectra acquired from dataset SN 26 and construct
the PFD curves obtained using both approaches. The good agreement between our
$E$-value based PFD (Sorić, 1989) and the target-decoy based PFD, displayed in
panel C of Figure 2, is expected because, as shown in panel D of Figure 1, we
have already found that the reported PFD and the observed PFD (computed by
using a known protein mixture) are nearly identical. The disagreement between
the $E$-value based PFD and the target-decoy based PFD using protocols of
(Spirin et al., 2011) seems to indicate that the moderate uncertainty in DPV
can influence the accuracy of the overall PFD estimate in a substantial
manner.

For RAId_DbS, the agreement between our $E$-value based PFD and the target-
decoy based PFD is further tested using more spectral datasets (SN16-SN25),
see supplementary Figure S6. In addition to RAId score, RAId_aPS allows other
scoring functions: XCorr, Hyperscore, and Kscore. For completeness, we plot
their corresponding protein $P$-value accuracy assessments in supplementary
Figures S7-S9; we also present the agreement tests between their $E$-value
based PFDs and the target-decoy based PFDs in supplementary Figures S10-S12.

###  7.3 Comparison with other methods

Figure 3:  Retrieval results of various methods based on their stated PFD
values. Because we can only ensure the accuracy of type-I error of the
proposed method, this figure only illustrates how our retrieval fares at the
stated value when compared to other methods. See the text of section 7.3 for
more details. To avoid clutter, results from using samples of Escherichia coli
whole cell lysyte, SN 27-30, are displayed in three panels (A, B, and C).
Within each panel, the results from RAId_DbS are always shown as a reference
curve. Similarly, results from using samples of Homo sapiens heart cells, SN
31-63, are also displayed in three panels (D, E, and F). The iProphet
(Shteynberg et al., 2011) results in panel E were downloaded from PeptideAtlas
instead of being computed.

The previous two subsections focus on the accuracy of type-I error control.
Although it is possible to accurately control type-I error for some protein ID
methods, this seems not the central focus of all protein ID methods. Many
protein ID methods prefer to use the decoy database search results to
pragmatically provide statistical significances for retrieval results from the
target (organismal) database. When this approach is used, the retrieval
results are displayed in terms of a parametric PFD plot: the parameter is some
kind of significance score used to prioritize the identifications, the
abscissa shows the proportion of false discovery and the ordinate displays the
number of identifications found in the target database. In general, a large
number of target identifications at a small PFD value indicates a good
retrieval, provided that the number of decoy identifications accurately
reflects the number of false identifications in the target database. However,
one should note that the fulfillment of the aforementioned condition requires
accurate type-I error control. Investigating and improving the statistical
accuracy of type-I error control of existing protein ID methods is beyond the
scope of this paper and we believe that it is best done by developers of
individual protein ID software.

To examine how our method compares with others under the pragmatic target-
decoy approach, we analyze two large datasets from E. Coli (SN 27-30) and Homo
Sapiens (SN 31-63) using a variety of protein ID software along with a number
of scoring functions. The list of software is given below (with both software
version and scoring functions, if given, shown inside a pair of parentheses):
RAId_DbS (v. Jan.12.2014; RAId), RAId_aPS (v. Jan.12.2014; XCorr, Kscore,
Hyperscore), Mascot (v. 2.4.0, http://www.matrixscience.com/help.html), and
X!Tandem (v. 2013.06.15; Hyperscore). The peptide identification software
SEQUEST (Eng et al., 1994) (v. 28) is only used in conjunction with other
post-processing protein ID software. We list below the post-processing
software used (with software version, peptide ID software, and peptide scoring
functions, if given, shown inside a pair of parentheses): iProphet (v. TPP
4.5; X!Tandem; Kscore), Proteome Discoverer (v. 1.3,
http://www.thermofisher.com
/en/home.html; SEQUEST, Mascot), and Scaffold Q+/Q+S (v. 4.0,
http://www.proteomesoftware.com; SEQUEST, Mascot). The results are displayed
in different panels of Figure 3. Before delving into the details of the
results, we first provide the information relevant to the generation of the
results.

In terms of peptide identification, RAId_DbS, RAId_aPS, Mascot, SEQUEST, and
X!Tandem used the same parameters: for Escherichia coli whole cell lysate, SN
27-30, the precursor ion mass error tolerance is $\pm$ 0.033 Da., the product
ion mass error tolerance is $\pm$ 0.033 Da., and up to 5 missed cleavages are
allowed; for Homo sapiens heart cells, SN 31-63, the precursor ion mass error
tolerance is $\pm$ 1.4 Da., the product ion mass error tolerance is $\pm$ 0.4
Da., and up to 2 missed cleavages are allowed.

Both X!Tandem and Mascot have built-in protein ID capability, and the target-
decoy approach was directly applied to estimate the protein level PFD. The
peptide ID outputs from SEQUEST and Mascot were also further analyzed using
Proteome Discoverer for protein identification and the target-decoy approach
was applied to estimate PFD. For iProphet, we did not compute the PFD but
downloaded the results for data group 4 from PeptideAtlas. Peptide
identification in this case was done using X!Tandem (v. 2009.10.01; Kscore).

Whenever the decoy peptide search results are available, Scaffold computes the
PFDs using the target-decoy approach; otherwise, it computes the PFDs using a
probabilistic method. In Figure 3 three Scaffold PFD curves are displayed, two
of which (shown in triangles and circles) are from target-decoy approaches.
The protein PFDs under Scaffold were computed by fixing the peptide threshold
at 20% PFD with a minimum of 1 evidence peptide per protein. We observed that
changing the peptide threshold to lower values had a small effect on the
number of proteins identified. We thus used the minimum of 1 peptide per
protein to maintain consistency across all methods. For RAId_DbS and RAId_aPS,
the PFD estimates do not require user-added target-decoy methods. RAId_DbS and
RAId_aPS compute the PFDs using the Sorić formula (Sorić, 1989).

Examinations of different panels of Figure 3 indicates that the retrieval
efficacy of the proposed method (shown in RAId_DbS and RAId_aPS PFD curves) is
comparable with existing protein ID methods, even though only at the stated
values. However, it should be noted that the proposed method does have a few
advantages. First, it reports accurate protein $P$-values, providing accurate
type-I error control. Second, the PFD curves obtained using this method show
stability across different mass resolution requirement and data sets, while
some methods seem to exhibit fluctuations of notable amplitudes.

##  8 Discussion

Our investigation indicates that it is possible to achieve faithful protein
$P$-value assignment, hence accurate type-I error control, in protein
identifications. Since our approach is founded on a derived mathematical
formula that requires accurate peptide $E$-values as input, it is evident that
accurate protein $P$-values require accurate statistical significance at the
peptide identification level.

The discrepancy between the computed protein $P$-value and the PFD results in
our implementation of the method of Spirin et al. (2011) is interesting. Based
on the results in Fig. 2, the peptide $P$-values are reasonably accurate
albeit exhibiting slightly larger fluctuations than the results from RAId_DbS.
In addition to the possibility of accumulating uncertainty of peptides’
$P$-values, the other possibility is that the iterative procedure to choose
the combination yielding the most significant Z-score may skew the $P$-values
towards the significant side. Investigation of the origin of the PFD and
$P$-value discrepancy when using the method of Spirin et al. (2011), however,
is beyond the scope of the current study and might be most appropriately done
by the authors of (Spirin et al., 2011).

As explained earlier, we allow more than one candidate peptide per spectrum to
accommodate the possibility of peptide co-elution. However, readers may ask
why do we choose to use DPVs for lower-ranking peptides per spectrum instead
of using ordered statistics. The reason is that in this context using ordered
statistics beyond the first is not meaningful: the $n$th ordered statistics
assumes that for a given query spectrum the best $n-1$ scored peptides are
spurious while the rank-$n$ peptide is the underlying peptide whose
fragmentation yields the query spectrum. This contradicts the general idea of
using a scoring function: among candidate peptides of a query spectrum, the
better a peptide scores the more likely it is the underlying peptide. On the
other hand, when using the DPV for the rank-$n$ peptide, we are essentially
assuming that the top $n-1$ candidate peptides of the query spectrum are co-
eluted underlying peptides and are not considered to be spurious.

The protein identification method proposed in this paper illustrates the
possibility of accurate type-I error control, providing a theoretically sound
significance assignment method that is also pragmatically simpler than the
target-decoy approach. This is particularly important since the number of
identified proteins versus PFDs provides trustworthy retrieval results only if
the reported PFDs truly reflects the proportion of false discoveries.
Evidently, to achieve accurate type-I error control is a task best done by
developers of individual software. Only when this is accomplished can a true
retrieval comparison among different methods be done.

Since we did not focus on type-II error, there is definite room for
improvement in terms of retrieval efficacy. We note that the information of
negatives (segments of a candidate protein not covered by the protein’s
evidence peptides) is not used. We also believe that, in principle, scoring
functions for peptide identification can also be improved to better separate
true underlying peptides from false positives. Currently, we are using a flat
peptide weight (by the number of proteins a peptide covers). It is perceivable
that more sophisticated weighting may be useful in better separating true
positive proteins from false positives. It is our plan to investigate these
avenues of improvement in the near future.

## Acknowledgement

#### Funding

This work was supported by the Intramural Research Program of the National
Library of Medicine at the National Institutes of Health.

## References

  * Alves and Yu (2008) Alves, G. and Yu, Y. K. (2008).  Statistical Characterization of a 1D Random Potential Problem - with applications in score statistics of MS-based peptide sequencing.  Physica A, 387(26), 6538–6544.
  * Alves and Yu (2011) Alves, G. and Yu, Y. K. (2011).  Combining independent, weighted P-values: achieving computational stability by a systematic expansion with controllable accuracy.  PLoS ONE, 6(8), e22647.
  * Alves et al. (2007) Alves, G., Ogurtsov, A. Y., and Yu, Y. K. (2007).  RAId_DbS: peptide identification using database searches with realistic statistics.  Biol. Direct, 2, 25.
  * Alves et al. (2008a) Alves, G., Ogurtsov, A. Y., Kwok, S., Wu, W. W., Wang, G., Shen, R. F., and Yu, Y. K. (2008a).  Detection of co-eluted peptides using database search methods.  Biol. Direct, 3, 27.
  * Alves et al. (2008b) Alves, G., Wu, W. W., Wang, G., Shen, R. F., and Yu, Y. K. (2008b).  Enhancing peptide identification confidence by combining search methods.  J. Proteome Res., 7(8), 3102–3113.
  * Alves et al. (2010) Alves, G., Ogurtsov, A. Y., and Yu, Y. K. (2010).  RAId_aPS: MS/MS analysis with multiple scoring functions and spectrum-specific statistics.  PLoS ONE, 5(11), e15438.
  * Bahrucha-Reid (1960) Bahrucha-Reid, A. (1960).  Elements of the Theory of Markov Processes and their Applications. McGraw-Hill.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.  Journal of the Royal Statistical Society. Series B (Methodological), 57(1), 289–300.
  * Elias and Gygi (2007) Elias, J. E. and Gygi, S. P. (2007).  Target-decoy search strategy for increased confidence in large-scale protein identifications by mass spectrometry.  Nat. Methods, 4(3), 207–214.
  * Eng et al. (1994) Eng, J. K., McCormack, A. L., and III, J. R. Y. (1994).  An approach to correlate tandem mass spectral data of peptides with amino acid sequences in a protein database.  Journal of the American Society for Mass Spectrometry, 5(11), 976 – 989.
  * Feng et al. (2007) Feng, J., Naiman, D. Q., and Cooper, B. (2007).  Probability-based pattern recognition and statistical framework for randomization: modeling tandem mass spectrum/peptide sequence false match frequencies.  Bioinformatics, 23(17), 2210–2217.
  * Fenyo and Beavis (2003) Fenyo, D. and Beavis, R. C. (2003).  A method for assessing the statistical significance of mass spectrometry-based protein identifications using general scoring schemes.  Anal. Chem., 75(4), 768–774.
  * Fenyo et al. (2010) Fenyo, D., Eriksson, J., and Beavis, R. (2010).  Mass spectrometric protein identification using the global proteome machine.  Methods Mol. Biol., 673, 189–202.
  * Fisher (1932) Fisher, R. A. (1932).  Statistical Methods for Research Workers, vol. II.  Oliver and Boyd, Edinburgh.
  * Good (1955) Good, I. J. (1955).  On the weighted combination of significance tests.  Journal of the Royal Statistical Society Series B (Methodological), 17, 264–265.
  * Gupta et al. (2011) Gupta, N., Bandeira, N., Keich, U., and Pevzner, P. A. (2011).  Target-decoy approach and false discovery rate: when things may go wrong.  J. Am. Soc. Mass Spectrom., 22(7), 1111–1120.
  * Huang et al. (2012) Huang, T., Wang, J., Yu, W., and He, Z. (2012).  Protein inference: a review.  Briefings in Bioinformatics.
  * Li and Radivojac (2012) Li, Y. F. and Radivojac, P. (2012).  Computational approaches to protein inference in shotgun proteomics.  BMC Bioinformatics, 13 Suppl 16, S4.
  * Mathai (1983) Mathai, A. (1983).  On linear combinations of independent exponential variables.  Communications in Statistics - Theory and Methods, 12(6), 625–632.
  * McHugh and Arthur (2008) McHugh, L. and Arthur, J. W. (2008).  Computational methods for protein identification from mass spectrometry data.  PLoS Comput. Biol., 4(2), e12.
  * Nesvizhskii et al. (2003) Nesvizhskii, A. I., Keller, A., Kolker, E., and Aebersold, R. (2003).  A statistical model for identifying proteins by tandem mass spectrometry.  Anal. Chem., 75(17), 4646–4658.
  * Noble and MacCoss (2012) Noble, W. S. and MacCoss, M. J. (2012).  Computational and statistical analysis of protein mass spectrometry data.  PLoS Comput. Biol., 8(1), e1002296.
  * Park et al. (2008) Park, C. Y., Klammer, A. A., Kall, L., MacCoss, M. J., and Noble, W. S. (2008).  Rapid and accurate peptide identification from tandem mass spectra.  J. Proteome Res., 7(7), 3022–3027.
  * Price et al. (2007) Price, T. S., Lucitt, M. B., Wu, W., Austin, D. J., Pizarro, A., Yocum, A. K., Blair, I. A., FitzGerald, G. A., and Grosser, T. (2007).  EBP, a program for protein identification using multiple tandem mass spectrometry datasets.  Mol. Cell Proteomics, 6(3), 527–536.
  * Robinson and Robinson (1991) Robinson, A. B. and Robinson, L. R. (1991).  Distribution of glutamine and asparagine residues and their near neighbors in peptides and proteins.  Proc. Natl. Acad. Sci. USA, 88, 8880–8884.
  * Searle (2010) Searle, B. C. (2010).  Scaffold: a bioinformatic tool for validating MS/MS-based proteomic studies.  Proteomics, 10(6), 1265–1269.
  * Segal (2008) Segal, M. R. (2008).  On E-values for tandem MS scoring schemes.  Bioinformatics, 24(14), 1652–1653.
  * Serang and Noble (2012) Serang, O. and Noble, W. (2012).  A review of statistical methods for protein identification using tandem mass spectrometry.  Stat Interface, 5(1), 3–20.
  * Serang et al. (2013) Serang, O., Paulo, J., Steen, H., and Steen, J. A. (2013).  A non-parametric cutout index for robust evaluation of identified proteins.  Mol. Cell Proteomics, 12(3), 807–812.
  * Shteynberg et al. (2011) Shteynberg, D., Deutsch, E. W., Lam, H., Eng, J. K., Sun, Z., Tasman, N., Mendoza, L., Moritz, R. L., Aebersold, R., and Nesvizhskii, A. I. (2011).  iProphet: multi-level integrative analysis of shotgun proteomic data improves peptide and protein identification rates and error estimates.  Mol. Cell Proteomics, 10(12), M111.007690.
  * Sorić (1989) Sorić, B. (1989).  Statistical “discoveries” and effect-size estimation.  Journal of the American Statistical Association, 84(406), 608–610.
  * Spirin et al. (2011) Spirin, V., Shpunt, A., Seebacher, J., Gentzel, M., Shevchenko, A., Gygi, S., and Sunyaev, S. (2011).  Assigning spectrum-specific P-values to protein identifications by mass spectrometry.  Bioinformatics, 27(8), 1128–1134.
  * Wang et al. (2009) Wang, G., Wu, W. W., Zhang, Z., Masilamani, S., and Shen, R. F. (2009).  Decoy methods for assessing false positives and false discovery rates in shotgun proteomics.  Anal. Chem., 81(1), 146–159.
  * Whitlock (2005) Whitlock, M. C. (2005).  Combining probability from independent tests: the weighted Z-method is superior to Fisher’s approach.  J. Evol. Biol., 18(5), 1368–1373.
  * Yang et al. (2004) Yang, X., Dondeti, V., Dezube, R., Maynard, D. M., Geer, L. Y., Epstein, J., Chen, X., Markey, S. P., and Kowalak, J. A. (2004).  DBParser: web-based software for shotgun proteomic data analyses.  J. Proteome Res., 3(5), 1002–1008.
  * Yu et al. (2006) Yu, Y. K., Gertz, E. M., Agarwala, R., Schaffer, A. A., and Altschul, S. F. (2006).  Retrieval accuracy, statistical significance and compositional similarity in protein sequence database searches.  Nucleic Acids Res., 34(20), 5966–5973.
  * Zhang and Chait (2000) Zhang, W. and Chait, B. T. (2000).  ProFound: an expert system for protein identification using mass spectrometric peptide mapping information.  Anal. Chem., 72(11), 2482–2489.
  * Zhang et al. (2013) Zhang, Y., Fonslow, B. R., Shan, B., Baek, M. C., and Yates, J. R. (2013).  Protein analysis by shotgun/bottom-up proteomics.  Chem. Rev., 113(4), 2343–2394.

## Supplementary Information

### Example application of formula (3) of the main text

To illustrate the use of the formula used to compute a protein $P$-value, let
us consider the following toy example. Let protein $\Pi$ have six peptide
evidences $\pi_{1},~{}\pi_{2},\ldots,\pi_{6}$ that falls into three groups
$\\{\pi_{1}\\}$, $\\{\pi_{2},\pi_{3}\\}$, and $\\{\pi_{4},\pi_{5},\pi_{6}\\}$,
respectively with weights $1$, $1/3$, $1/2$. This means that peptide $\pi_{1}$
is a subsequence of protein $\Pi$ only, peptides $\pi_{2}$ and $\pi_{3}$ are
subsequences of three proteins ($\Pi$ and two others), and peptides $\pi_{4}$,
$\pi_{5}$, and $\pi_{6}$ are subsequences of two proteins ($\Pi$ and another
one). Let the peptide $E$-values of these six evidence peptides be
$e_{1},~{}e_{2},\ldots,~{}e_{6}$, all less than one. Under our approximating
database $P$-value by $E$-value, this means that the evidence peptides have
their respective database $P$-values $e_{1},~{}e_{2},\ldots,~{}e_{6}$. From
the information above we know that $m=3$, $n_{1}=1$, $n_{2}=2$, $n_{3}=3$,
$w_{1}=1$, $w_{2}=1/3$, $w_{3}=1/2$, $r_{1}=1/w_{1}=1$, $r_{2}=1/w_{2}=3$, and
$r_{3}=1/w_{3}=2$.

To use formula (3) in the manuscript to compute a protein $P$-value, we first
need the quantity $\tau$ given by eq. (1) in the main text. In the current
case– $m=3$, $n_{1}=1$, $n_{2}=2$, $n_{3}=3$ and with peptides database
$P$-values $e_{1},~{}e_{2},\ldots,~{}e_{6}$ – the quantity $\tau$ can be
written as

$\tau=\prod_{k=1}^{3}\left[\prod_{j=1}^{n_{k}}e_{k;j}\right]^{w_{k}}=e_{1}^{w_{1}}e_{2}^{w_{2}}e_{3}^{w_{2}}e_{4}^{w_{3}}e_{5}^{w_{3}}e_{6}^{w_{3}}.$

And the protein $P$-value is given by

$F(\tau)=\left[\prod_{l=1}^{3}r_{l}^{n_{l}}\right]\sum_{k=1}^{3}\,\sum_{\mathcal{G}(k)}\;\biggr{\\{}\frac{1}{r_{k}^{g_{k}+1}}H(-r_{k}\ln\tau\,,\,g_{k})\times\biggr{(}\prod_{j=1,j\neq
k}^{3}\frac{(n_{j}-1+g_{j})!}{(n_{j}-1)!g_{j}!}\frac{(-1)^{g_{j}}}{(r_{j}-r_{k})^{n_{j}+g_{j}}}\biggr{)}\biggr{\\}},$
(6)

where $H(-r_{k}\ln\tau\,,\,g_{k})$ is given by

$H(-r_{k}\ln\tau\,,\,g_{k})=\exp(r_{k}\ln\tau)\sum_{k=0}^{g_{k}}\frac{(-r_{k}\ln\tau)^{k}}{k!}.$

Remember that $\sum_{\mathcal{G}(k)}$ enumerates each set of nonnegative
integers $\\{g_{1},g_{2},\ldots,g_{m}\\}$ that satisfies the $k$-dependent
constraint $\sum_{i=1}^{m}g_{i}=n_{k}-1$, it is thus possible to replace the
$\sum_{\mathcal{G}(k)}$ in eq. (6) by an $m$-dimensional summation with an
explicit constraint. Specifically, we can rewrite the sum over set as

$\displaystyle F(\tau)$ $\displaystyle=$
$\displaystyle\left[\prod_{l=1}^{3}r_{l}^{n_{l}}\right]\sum_{k=1}^{3}\,\sum_{g_{1}=0}^{n_{k}-1}\sum_{g_{2}=0}^{n_{k}-1}\sum_{g_{3}=0}^{n_{k}-1}\delta_{g_{1}+g_{2}+g_{3},n_{k}-1}\;\biggr{\\{}\frac{1}{r_{k}^{g_{k}+1}}H(-r_{k}\ln\tau\,,\,g_{k})\times$
(7) $\displaystyle\times\biggr{(}\prod_{j=1,j\neq
k}^{3}\frac{(n_{j}-1+g_{j})!}{(n_{j}-1)!g_{j}!}\frac{(-1)^{g_{j}}}{(r_{j}-r_{k})^{n_{j}+g_{j}}}\biggr{)}\biggr{\\}}.$

where the Kronecker delta function $\delta_{n,n^{\prime}}$ takes value one if
$n=n^{\prime}$ but zero otherwise.

The first product on the right hand side of eq. (7) is equal to

$\left[\prod_{l=1}^{3}r_{l}^{n_{l}}\right]=(r_{1})^{1}\cdot(r_{2})^{2}\cdot(r_{3})^{3}=1^{1}\cdot
3^{2}\cdot 2^{3}\,.$

After this overall factor is obtained, the main task is to evaluate the
summation over $k$, which ranges from $k=1$ to $k=3$. For each $k$, we are
only interested in the non-negative integral $g_{i}$s that satisfy the
$k$-dependent constraint $\sum_{i=1}^{3}g_{i}=n_{k}-1$. When $k=1$ we have
$n_{1}=1$ and the constrained summation

$\sum_{g_{1}=0}^{0}\sum_{g_{2}=0}^{0}\sum_{g_{3}=0}^{0}\delta_{g_{1}+g_{2}+g_{3},0}$

only allows one valid {$g_{1},g_{2},g_{3}$} set, namely, {0,0,0}. For $k=2$ we
have $n_{2}=2$ and the constrained summation

$\sum_{g_{1}=0}^{1}\sum_{g_{2}=0}^{1}\sum_{g_{3}=0}^{1}\delta_{g_{1}+g_{2}+g_{3},1}$

allows three sets of valid {$g_{1},g_{2},g_{3}$}, namely, {1,0,0}, {0,1,0} and
{0,0,1}. For $k=3$ we have $n_{3}=3$ and the constrained summation

$\sum_{g_{1}=0}^{2}\sum_{g_{2}=0}^{2}\sum_{g_{3}=0}^{2}\delta_{g_{1}+g_{2}+g_{3},2}$

allows six sets of valid {$g_{1},g_{2},g_{3}$}, namely, {1,1,0}, {1,0,1},
{0,1,1}, {2,0,0}, {0,2,0} and {0,0,2}. Each valid set of
$\\{g_{1},g_{2},g_{3}\\}$ must be substituted into the summand (inside the
pair of curly braces) of eq. (7) to yield its respective contribution for the
$P$-value.

### Supplementary Tables

Table S1:  MS/MS data group 1. The corresponding spectral datasets are produced by using Sigma49 (a protein standard mixture composed of 49 known human proteins), and are downloaded from the National Center for Biotechnology Information (Peptidome database) at ftp://ftp.ncbi.nih.gov/pub/peptidome/studies/PSEnnn/PSE108ftp://ftp.ncbi.nih.gov/pub/peptidome/studies/PSEnnn/PSE108. The column heading SN represents sample index, the abbreviation CGL stands for chromatography gradient length, and the column heading $n_{s}$ stands for the number of spectra acquired. The rest of the column headings are self-explanatory.  SN | Sample Load | Instrument | CGL(minutes) | $n_{s}$ | File Name
---|---|---|---|---|---
1 | 5 fmol | LTQ Orbitrap | 45 | 1,531 | PSM1027_07FEB15_ABRF_FT_5a.mzXML
2 | 5 fmol | LTQ Orbitrap | 45 | 1,902 | PSM1028_07FEB15_ABRF_FT_5b.mzXML
3 | 5 fmol | LTQ Orbitrap | 45 | 2,014 | PSM1029_07FEB15_ABRF_FT_5c.mzXML
4 | 10 fmol | LTQ Orbitrap | 45 | 2,026 | PSM1027_07FEB15_ABRF_FT_10a.mzXM
5 | 10 fmol | LTQ Orbitrap | 45 | 2,125 | PSM1028_07FEB15_ABRF_FT_10b.mzXML
6 | 10 fmol | LTQ Orbitrap | 45 | 2,253 | PSM1029_07FEB15_ABRF_FT_10c.mzXML
7 | 25 fmol | LTQ Orbitrap | 45 | 2,772 | PSM1027_07FEB15_ABRF_FT_25a.mzXML
8 | 25 fmol | LTQ Orbitrap | 45 | 2,669 | PSM1028_07FEB15_ABRF_FT_25b.mzXML
9 | 25 fmol | LTQ Orbitrap | 45 | 2,504 | PSM1029_07FEB15_ABRF_FT_25c.mzXML
10 | 50 fmol | LTQ Orbitrap | 45 | 3,259 | PSM1027_07FEB15_ABRF_FT_50a.mzXML
11 | 50 fmol | LTQ Orbitrap | 45 | 3,406 | PSM1028_07FEB15_ABRF_FT_50b.mzXML
12 | 50 fmol | LTQ Orbitrap | 45 | 2,993 | PSM1029_07FEB15_ABRF_FT_50c.mzXML
13 | 100 fmol | LTQ Orbitrap | 45 | 3,629 | PSM1027_07FEB15_ABRF_FT_100a.mzXML
14 | 100 fmol | LTQ Orbitrap | 45 | 3,622 | PSM1028_07FEB15_ABRF_FT_100b.mzXML
15 | 100 fmol | LTQ Orbitrap | 45 | 3,592 | PSM1029_07FEB15_ABRF_FT_100c.mzXML
Table S2:  MS/MS data group 2. The corresponding spectral datasets are produced by using a complex protein mixture of Escherichia coli K-12 whole cell lysate, and are downloaded from the Pacific Northwest National Laboratory at http://omics.pnl.gov/http://omics.pnl.gov/.  SN | Instrument | CGL(minutes) | $n_{s}$ | File Name
---|---|---|---|---
16 | LTQ Orbitrap | 100 | 18,573 | Ecoli432_R1-rr_18Dec09_Falcon_09-09-14.mzXML
17 | LTQ Orbitrap | 100 | 18,585 | Ecoli432_R2_7Dec09_Falcon_09-09-15.mzXML
18 | LTQ Orbitrap | 100 | 18,669 | Ecoli432_R3_7Dec09_Falcon_09-09-16.mzXML
19 | LTQ Orbitrap | 100 | 18,585 | Ecoli432_R4_15Dec09_Falcon_09-09-16.mzXML
20 | LTQ Orbitrap | 100 | 18,650 | Ecoli433_R1_7Dec09_Falcon_09-09-14.mzXML
21 | LTQ Orbitrap | 100 | 18,763 | Ecoli433_R2_7Dec09_Falcon_09-09-15.mzXML
22 | LTQ Orbitrap | 100 | 18,770 | Ecoli433_R4_13Dec09_Falcon_09-09-16.mzXML
23 | LTQ Orbitrap | 100 | 18,488 | Ecoli434_R1_7Dec09_Falcon_09-09-14.mzXML
24 | LTQ Orbitrap | 100 | 18,923 | Ecoli434_R2_7Dec09_Falcon_09-09-15.mzXML
25 | LTQ Orbitrap | 100 | 19,010 | Ecoli434_R3_7Dec09_Falcon_09-09-16.mzXML
26 | LTQ Orbitrap | 100 | 18,737 | Ecoli434_R4_13Dec09_Falcon_09-09-16.mzXML
Table S3:  MS/MS data group 3. The corresponding spectral datasets are produced by using a complex protein mixture of Escherichia coli whole cell lysate prepared in house.  SN | Instrument | CGL(minutes) | $n_{s}$ | File Name
---|---|---|---|---
27 | Orbitrap Elite | 90 | 24,280 | E_L_2.mzML
28 | Orbitrap Elite | 90 | 22,435 | E_M_2.mzML
29 | Orbitrap Elite | 90 | 23,875 | E_H_2.mzML
30 | Orbitrap Elite | 90 | 18,573 | E_S_2.mzML
Table S4:  MS/MS data group 4. The corresponding spectral datasets are produced by using SDS-PAGE protein fractionation extraction of human lung cells, and are downloaded from PeptideAtlas database at ftp://ftp.peptideatlas.org/pub/PeptideAtlas/Repository/PAe001771ftp://ftp.peptideatlas.org/pub/PeptideAtlas/Repository/PAe001771.  SN | Instrument | CGL(minutes) | $n_{s}$ | File Name
---|---|---|---|---
31-63 | LTQ Orbitrap | 120 | 340,861 | Roche_human_lung_001.mzML - Roche_human_lung_033.mzML

### Supplementary Figures

Figure S1:  Accuracy assessment of the protein P-value reported by the
proposed method using data group 1 and RAId_DbS’s peptide $E$-values. All
spectra in data group 1 are used to search the Homo sapiens database with
precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion mass tolerance $\pm$
0.8 Da. and up-to-2 missed cleavage sites allowed. The accuracy of the
reported protein P-value is evaluated by plotting it versus the observed
protein P-value, see main text for details. The closer the above curves are to
the $y=x$ line the more accurate are the reported protein P-values. The two
dash lines, $y=2x$ and $y=0.5x$, are provided as visual guides regarding how
close/off the reported protein P-values are to the $y=x$ line.

Figure S2:  Accuracy assessment of the protein P-value reported by the
proposed method using data group 2 and RAId_DbS’s peptide $E$-values. All
spectra in data group 2 are used to search the Escherichia coli database with
precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion mass tolerance $\pm$
0.033 Da. and up-to-2 missed cleavage sites allowed. The accuracy of the
reported protein P-value is evaluated by plotting it versus the observed
protein P-value, see main text for details. The closer the above curves are to
the $y=x$ line the more accurate are the reported protein P-values. The two
dash lines, $y=2x$ and $y=0.5x$, are provided as visual guides regarding how
close/off the reported protein P-values are to the $y=x$ line.

Figure S3:  Accuracy assessment of the reported protein P-value by the
proposed method when the mass tolerances are very small. The protein
$P$-values are obtained by using RAId_DbS’s peptide $E$-values. All spectra in
data group 2 are used to search the Escherichia coli database with precursor-
ion mass tolerance $\pm$ 0.0033 Da., product-ion mass tolerance $\pm$ 0.0033
Da. and up-to-2 missed cleavage sites allowed. The accuracy of the reported
protein P-value is evaluated by plotting it versus the observed protein
P-value, see main text for details. The closer the above curves are to the
$y=x$ line the more accurate are the reported protein P-values. The two dash
lines, $y=2x$ and $y=0.5x$, are provided as visual guides regarding how
close/off the reported protein P-values are to the $y=x$ line.

Figure S4: Accuracy assessment of the protein E-value reported by the proposed
method using data group 1 and RAId_DbS’s peptide $E$-values. All spectra in
data group 1 are used to search the Homo sapiens database with precursor-ion
mass tolerance $\pm$ 0.033 Da., product-ion mass tolerance $\pm$ 0.8 Da. and
up-to-2 missed cleavage sites allowed. The accuracy of the reported protein
P-value is evaluated by plotting it versus the observed protein P-value, see
main text for details. The closer the above curves are to the $y=x$ line the
more accurate are the reported protein E-values. The two dash lines, $y=2x$
and $y=0.5x$, are provided as visual guides regarding how close/off the
reported protein E-values are to the $y=x$ line.

Figure S5:  Accuracy assessment of reported protein PFD using data group 1 and
RAId_DbS’s peptide $E$-values. All spectra in data group 1 (MS/MS spectra from
SN 1-15) are used to search the Homo sapiens database with precursor-ion mass
tolerance $\pm$ 0.033 Da., product-ion mass tolerance $\pm$ 0.8 Da. and up-
to-2 missed cleavage sites allowed. The reported PFD is computed by using the
reported protein E-value in Sorić’s formula, while the observed PFD is
obtained from the ratio of the number of false discoveries (false positives)
to the total number of discoveries (true positives + false positives) at a
given E-value cutoff. The closer the observed and reported curves are to each
other the better.

Figure S6:  The agreement between the $E$-value based PFD and the target-decoy
based PFD when the peptide $E$-values are from RAId_DbS. All spectra in data
group 2 (MS/MS spectra SN16-26) are used to search the Escherichia coli
database with precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion mass
tolerance $\pm$ 0.033 Da. and up-to-2 missed cleavage sites allowed. The
$E$-value based PFD is computed by using the reported protein E-value in
Sorić’s formula, while the target-decoy based PFD is obtained from the ratio
of the number of identified decoy proteins to the total number of identified
proteins (target + decoy) for a given E-value cutoff. Within each panel, the
closer the two curves are to each other the better.

Figure S7:  Accuracy assessment of the protein P-value using data group 2 and
RAId_aPS’s peptide $E$-values when the selected scoring function is XCorr. All
spectra in data group 2 are used to search the Escherichia coli database with
precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion mass tolerance $\pm$
0.033 Da. and up-to-2 missed cleavage sites allowed. The accuracy of the
reported protein P-value is evaluated by plotting it versus the observed
protein P-value, see main text for details. The closer the above curves are to
the $y=x$ line the more accurate are the reported protein P-values. The two
dash lines, $y=2x$ and $y=0.5x$, are provided as visual guides regarding how
close/off the reported protein P-values are to the $y=x$ line.

Figure S8: Accuracy assessment of the protein P-value using data group 2 and
RAId_aPS’s peptide $E$-values when the selected scoring function is
Hyperscore. All spectra in data group 2 are used to search the Escherichia
coli database with precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion
mass tolerance $\pm$ 0.033 Da. and up-to-2 missed cleavage sites allowed. The
accuracy of the reported protein P-value is evaluated by plotting it versus
the observed protein P-value, see main text for details. The closer the above
curves are to the $y=x$ line the more accurate are the reported protein
P-values. The two dash lines, $y=2x$ and $y=0.5x$, are provided as visual
guides regarding how close/off the reported protein P-values are to the $y=x$
line.

Figure S9:  Accuracy assessment of the protein P-value using data group 2 and
RAId_aPS’s peptide $E$-values when the selected scoring function is Kscore.
All spectra in data group 2 are used to search the Escherichia coli database
with precursor-ion mass tolerance $\pm$ 0.033 Da., product-ion mass tolerance
$\pm$ 0.033 Da. and up-to-2 missed cleavage sites allowed. The accuracy of the
reported protein P-value is evaluated by plotting it versus the observed
protein P-value, see main text for details. The closer the above curves are to
the $y=x$ line the more accurate are the reported protein P-values. The two
dash lines, $y=2x$ and $y=0.5x$, are provided as visual guides regarding how
close/off the reported protein P-values are to the $y=x$ line.

Figure S10:  The agreement between the $E$-value based PFD and the target-
decoy based PFD when the peptide $E$-values are from RAId_aPS(XCorr). All
spectra in data group 2 (MS/MS spectra SN16-26) are used to search the
Escherichia coli database with precursor-ion mass tolerance $\pm$ 0.033 Da.,
product-ion mass tolerance $\pm$ 0.033 Da. and up-to-2 missed cleavage sites
allowed. The $E$-value based PFD is computed by using the reported protein
E-value in Sorić’s formula, while the target-decoy based PFD is obtained from
the ratio of the number of identified decoy proteins to the total number of
identified proteins (target + decoy) for a given E-value cutoff. Within each
panel, the closer the two curves are to each other the better.

Figure S11:  The agreement between the $E$-value based PFD and the target-
decoy based PFD when the peptide $E$-values are from RAId_aPS(Hyperscore). All
spectra in data group 2 (MS/MS spectra SN16-26) are used to search the
Escherichia coli database with precursor-ion mass tolerance $\pm$ 0.033 Da.,
product-ion mass tolerance $\pm$ 0.033 Da. and up-to-2 missed cleavage sites
allowed. The $E$-value based PFD is computed by using the reported protein
E-value in Sorić’s formula, while the target-decoy based PFD is obtained from
the ratio of the number of identified decoy proteins to the total number of
identified proteins (target + decoy) for a given E-value cutoff. Within each
panel, the closer the two curves are to each other the better.

Figure S12:  The agreement between the $E$-value based PFD and the target-
decoy based PFD when the peptide $E$-values are from RAId_aPS(Kscore). All
spectra in data group 2 (MS/MS spectra SN16-26) are used to search the
Escherichia coli database with precursor-ion mass tolerance $\pm$ 0.033 Da.,
product-ion mass tolerance $\pm$ 0.033 Da. and up-to-2 missed cleavage sites
allowed. The $E$-value based PFD is computed by using the reported protein
E-value in Sorić’s formula, while the target-decoy based PFD is obtained from
the ratio of the number of identified decoy proteins to the total number of
identified proteins (target + decoy) for a given E-value cutoff. Within each
panel, the closer the two curves are to each other the better.


11institutetext: German Cancer Research Center (DKFZ), Heidelberg, Division of
Medical Image Computing, Germany
22institutetext: Helmholtz Imaging, DKFZ, Heidelberg, Germany
33institutetext: Pattern Analysis and Learning Group, Department of Radiation
Oncology, Heidelberg University Hospital, Heidelberg, Germany
44institutetext: National Center for Tumor Diseases (NCT), NCT Heidelberg, A
partnership between DKFZ and University Medical Center Heidelberg
55institutetext: Medical Faculty Heidelberg, University of Heidelberg,
Heidelberg, Germany 66institutetext: Faculty of Mathematics and Computer
Science, Heidelberg University, Germany 66email: constantin.ulrich@dkfz-
heidelberg.de

# Mitigating False Predictions In Unreasonable Body Regions

Constantin Ulrich 114455 Catherine Knobloch 11 Julius C. Holzschuh 11 Tassilo
Wald 112266 Maximilian R. Rokuss 1166 Maximilian Zenk 1155 Maximilian Fischer
11 Michael Baumgartner 112266 Fabian Isensee 1122 Klaus H. Maier-Hein 1133

###### Abstract

Despite considerable strides in developing deep learning models for 3D medical
image segmentation, the challenge of effectively generalizing across diverse
image distributions persists. While domain generalization is acknowledged as
vital for robust application in clinical settings, the challenges stemming
from training with a limited Field of View (FOV) remain unaddressed. This
limitation leads to false predictions when applied to body regions beyond the
FOV of the training data. In response to this problem, we propose a novel loss
function that penalizes predictions in implausible body regions, applicable in
both single-dataset and multi-dataset training schemes. It is realized with a
Body Part Regression model that generates axial slice positional scores.
Through comprehensive evaluation using a test set featuring varying FOVs, our
approach demonstrates remarkable improvements in generalization capabilities.
It effectively mitigates false positive tumor predictions up to 85% and
significantly enhances overall segmentation performance.

###### Keywords:

Medical image segmentation generalization domain shift robustness foundation
model partially labeled datasets.

**footnotetext: Equal Contribution.

##  1 Introduction

Figure 1: a) shows an example of a training dataset with a limited Fiel of
View (FOV) caused by the desire to reduce radiation dose, scanning time and
costs. Typically, the final model is then applied to images from different
sources and studies with varying FOVs and tends to make anatomical implausible
predictions in unseen body parts. b) shows the FOV of the datasets of the
Medical Decathlon as well as the distribution of the location of the upper and
lower position of the target structure, which was determined using a Body Part
Regression model.

In modern healthcare, 3D medical image segmentation is crucial to support
clinicians in treatment planning, disease monitoring, and anatomical as well
as pathological structure analysis. Despite significant advancements of deep
learning models for medical image segmentation [10, 19, 9, 22, 11], a
persistent challenge remains: a robust generalization to image distribution
shifts. Many existing methodologies have attempted to address this by focusing
on improving generalization between clinics, scanner types, imaging protocols,
and demographics [4, 7, 13, 28]. However, a fundamental aspect often
overlooked is the limited Field of View (FOV) inherent in most 3D medical
image datasets used for model training. This restricted FOV primarily arises
due to constraints imposed by factors such as the reduction of radiation dose,
acquisition time, and associated costs, especially in datasets restricted to
certain pathologies of interest. Even if the images of a public dataset come
from different sources with originally larger and more diverse FOVs, the
images are often cropped during the manual annotation and publishing process
to streamline image processing and minimize memory consumption. However,
during inference, this results in tangible scenarios in which trained models
are applied to images with a different FOV. This is a particularly severe
problem when large image cohorts need to be automatically delineated.
The disparity between the restricted FOV during training and larger and more
diverse FOVs encountered in real-world applications poses a major challenge to
the adaptability and generalizability of segmentation models. As illustrated
in fig. 1 a), models tend to make incorrect predictions in unseen regions not
represented in the model’s training data e.g. entail liver vessels in the
head. Such false predictions are completely irrational for humans and these
conspicuous errors not only diminish the measurable performance of these
models but also drastically undermine the trust that patients and clinicians
place in these methods.
In contrast to previous methods that addressed false predictions by an error
prone additional post-processing step [12, 8, 23], we introduce a novel
training approach designed to mitigate wrong predictions in body parts without
label information. A novel Region Loss directly penalizes predictions in
previously defined unreasonable body regions during model training. To
determine the corresponding body regions, we employ a self-supervised trained
Body Part Regression (BPR) model, which assigns each axial CT image slice a
specific value corresponding to its axial position, ensuring that anatomically
equivalent body parts are assigned similar values across patients and images
[23, 27].
We evaluate the proposed method in both, a standard single-dataset (SD)
training scenario as well as a multi-dataset (MD) training scenario, where
multiple datasets with various FOVs and heterogenous class annotations are
used concurrently [24]. Importantly, as we illustrate the potential of our
method for 3D image segmentation in CT images, its versatility and the self-
supervised training of the BPR model pave the way for straightforward
expansion to a range of other 3D image modalities and tasks, including
detection. Our results demonstrate that the proposed method effectively allows
generalization between varying FOVs and mitigates false positive predictions
in anatomically unreasonable body parts. The code and the models will be
published at: tba.

##  2 Method

In the following, we present our approach to prevent predictions in
unreasonable body parts using a novel Region Loss. The method is implemented
both in a MD and in a SD setting in conjuction with an unlabeled support
dataset to provide the network with varying FOVs.

###  2.1 Body part Regression Model

To translate each anatomical body part of a CT slice into a standardized
machine-interpretable position we employed the publicly available BPR model by
Schuhegger et al. which is an extension of the work of Yan et al. [23, 27].
While the training details can be found in the corresponding publications, the
final model $f_{\theta}$ with parameters $\theta$ has learned a mapping of an
axial slice $\hat{x}$ of a CT image $X$ to a slice score $s$ where the start
of the pelvis is mapped to 0 and the head to 100:
$f_{\theta}:\mathbb{R}^{2}\xrightarrow{}\mathbb{R}:f_{\theta}(\hat{x})=s$. The
BPR model was applied to all training images to assign each image voxel
$x_{n}$ in each slice to the corresponding encoded axial position $s_{n}$. We
approximated the BPR scores for slices below the pelvis using linear
extrapolation. Fig. 1 (b) shows the image FOV, as well as the foreground class
distributions of all CT datasets of the Medical Segmentation Decathlon [2]).
None of the datasets covers all body parts from the pelvis to the head and
some datasets exclusively cover a FOV limited to the extent of the target
structure.

###  2.2 Region Loss

For training the segmentation networks, we introduce a novel loss function to
penalize predictions in unreasonable body parts. The loss is based on the
weighted Binary Cross Entropy (wBCE) loss defined for the network prediction
$\hat{y}$ and the corresponding binary ground truth $y$ of an individual voxel
as follows:

$wBCE(\hat{y}_{n},y_{n},w_{n})=w_{n}[y_{n}\log(\hat{y}_{n})+(1-y_{n})\log(1-\hat{y}_{n})]=w_{n}\log(1-\hat{y}_{n})=l_{n},$
(1)

where $w_{n}$ is a weighting factor for each predicted voxel $\hat{y}_{n}$
with $n\in[1,N]$ where N is the number of voxels in an image. This factor is
used to determine in which body part the loss term penalizes the predictions
and in which body part predictions were considered to be valid. Within the
invalid region, the ground truth of a target class is expected to be zero,
which translates to $y_{n}=0$. In section 2.3.4, it is explained how the upper
and lower boundaries for the valid body regions $S^{max},S^{min}$ were
determined. Taking into account the uncertainties of the BPR mapping and
variations in the human anatomy, we decided to choose a smooth transition of
the voxel weighting between valid and invalid regions using a Gaussian
smoothing:

$w_{n}(s_{n},S^{min},S^{max})=\begin{cases}0&:S^{min}\leq s_{n}\leq S^{max}\\\
1-\exp(\frac{-(s_{n}-S^{min})^{2}}{2(\sigma^{min})^{2}})&:s_{n}<S^{min}\\\
1-\exp(\frac{-(s_{n}-S^{max})^{2}}{2(\sigma^{max})^{2}})&:s_{n}>S^{max},\end{cases}$
(2)

where $\sigma^{max}$ and $\sigma^{min}$ correspond to the standard deviation
of the distribution of all minimum and maximum slice scores of a target class
over all training images.
The final Region Loss only penalizes the largest k percent of the wBCE loss
values of all predictions [17] made in the invalid regions. Otherwise, many
small predicted probabilities would reduce the overall loss signal:

$RL=TopK(L)=\alpha\frac{1}{M}\sum_{n=1}^{M}l_{n},$ (3)

where $\alpha$ is a weighting factor, $M=\hat{N}\cdot\frac{k}{100}$,
$L=\\{l_{n}|\forall n\in[1,\hat{N}]\\}$, $\hat{N}$ is the number of voxels
within the defined invalid region, and $l_{n}$ is the $n$-largest wBCE voxel
value.

###  2.3 Experiments

####  2.3.1 Multi-dataset training

In this study, we employ the loss function described above within the recently
introduced MultiTalent framework for MD training [24]. In this MD setting, a
network undergoes training using multiple datasets, each containing a distinct
subset of annotated target structures. With the substantial expenses
associated with annotating 3D medical images, the scientific community has
recently intensified its focus on the challenge of learning from partially
labeled datasets [18]. MultiTalent comes with a dataset and target adaptive
loss function to calculate the loss solely for the classes that were annotated
in the corresponding dataset, while simultaneously predicting all classes from
all datasets [24]. Consequently, the network encounters FOVs where ground
truth annotations for certain classes are absent. For these classes, we
incorporate the Region Loss penalize predictions in unreasonable body regions.
The total loss can be formulated as follows:

$L^{(k)}=\sum_{c^{(k)}}L_{c}+\sum_{c^{(\neg k)}}RL_{c}$ (4)

where we add up the class and adaptive loss $L$ over all annotated classes
$c^{(k)}\in C^{(k)}$ of a dataset $k\in K$ and the Region Loss $RL$ over all
unlabeled classes $c^{(\neg k)}\in C^{(\neg k)}$, which corresponds to all
classes that are annotated in the remaining datasets of $K\setminus k$. $K$
denoted the collection of datasets used for training the MultiTalent network.
The MultiTalent network was trained on 13 datasets as the original work [1,
15, 6, 16, 14, 20, 3, 21, 8]. We restricted the evaluation to the classes of 5
Medical Segmentation Decathlon (MSD) CT datasets [2].

####  2.3.2 Single-dataset training

In a conventional training scenario, where a network is trained using a SD
containing a specific set of annotated target structures, the network does not
encounter FOVs lacking ground truth annotations. Consequently, we chose 120 CT
images from the AutoPET2 dataset [5] that do not contain annotations of the
target classes but possesses a distinct and larger FOV. During training, we
select these images with a probability of $p$. For all annotated images from
the SD, the default nnU-Net supervised loss was computed [10], while for
images lacking label information, the Region Loss was employed. In this work,
we tested the proposed method for 5 MSD CT datasets [2].

####  2.3.3 Implementation

All experiments conducted in this study utilize the established nnU-Net
framework [10]. The slice scores are obtained by employing the BPR model
developed by Schuhegger et al. prior to training the segmentation model to
mitigate computational overhead. A three-dimensional array, which has the same
shape as the raw data and contains the slice scores for each corresponding
voxel at the same position, is generated after applying the BPR model to the
image. All spatial preprocessing steps and augmentations such as rotation,
translation, and mirroring, are applied to both the training data as well as
the slice score arrays. During training, these arrays are utilized to
determine the anatomical region of the input patches in which the Region Loss
penalizes predictions. While many design decisions align with the default
suggestions provided by nnU-Net and MultiTalent, an overview of all altered
hyperparameters can be found in the appendix in table 3.

####  2.3.4 Defining valid regions

To find the best boundaries that delineate the anatomical regions where
predictions should be penalized by Region Loss, multiple MultiTalent networks
were trained with incrementally increasing valid regions. Initially, we chose
the mean value of the minimum and maximum slice scores from the distribution
corresponding to the lowest and highest ground truth label positions for each
class. Because the distributions’ positions and widths differ across classes,
we extended the valid region step-wise according to the standard deviation of
the minimum and maximum slice score distribution. The outcomes for all classes
within the MultiTalent dataset collection are shown in the appendix in table
3. Given the variations of the best boundary definition across classes based
on the cross-validation results, we selected the best boundaries for each
class for all subsequent experiments.

####  2.3.5 Baselines

For both the MD and SD training settings, we employ the default training
approach alongside two additional postprocessing schemes as baselines. The
first postprocessing baseline utilizes the BPR model to retrospectively crop
all predictions outside of the valid regions (BPR postprocessing).
Consequently, the BPR model must first be applied to all test images to create
a mask for the valid regions. The boundaries for this postprocessing cropping
were determined through optimization on the validation predictions of the MD
training and the SD training. In the second postprocessing baseline, all
connected components except for the largest component LC were removed (LC
postprocessing). This approach is particularly applicable in organ
segmentation tasks where only one target instance is typically expected,
making it a well-established baseline. Our method stands as the first approach
that addresses the FOV issue directly during the training not relying on an
additional postprocessing step.

##  3 Results

Figure 2: The proposed Region Loss achieves always the best (red) or second
best (blue) mean Dice. Additionally, the Region Loss improves the boundaries
of the lower two boxplot quartiles compared to default MultiTalent or default
nnU-Net. This indicates, that the Region Loss succesfully improves inferior
cases.

To evaluate the efficacy of the proposed method, we used all 146 images from
the official test and validation set of the TotalsegmentatorV2 dataset [25].
This dataset includes comprehensive field-of-views (FOVs) and corresponding
organ annotations, akin to the annotations found in the MSD datasets. In fig.
2 on the left, we present the results of the MD training setting. Notably, the
proposed Region Loss consistently resulted in a better performance compared to
the default MultiTalent baseline across all datasets. The mean Dice score has
not only improved, but the narrower shape of the performance distributions,
visualized by the boxplots, also suggests an enhanced consistency in the
results with fewer negative outliers. Adding a BPR postprocessing step
enhanced pancreas and in particular liver segmentation. The boundaries for
cropping the predictions in this additional postprocessing step were optimized
on the 5-fold cross-validation results. When the unknown inference data
distribution resembles the training data, this postprocessing baseline serves
as an upper limit for our method. Conversely, applying BPR postprocessing
adversely impacted spleen segmentation, resulting in the removal of some true
positive predictions. In general, this postprocessing heavily relies on the
accuracy of the BPR model and the selection of optimal boundaries, which is
challenging without access to the inference data. Removing all but the largest
component only resulted in a notable positive impact on spleen segmentation.
In the SD results presented in fig. 2 on the right, the overall trend remains
consistent, and the proposed method shows a slight improvement for all organs
compared to default nnU-Net training. It is noteworthy that the MD results are
significantly better than the SD results. While the initial MultiTalent
publication hinted at a modest advantage over SD training within in-
distribution settings, this out-of-distribution scenario underscores
MultiTalent’s superior generalization capabilities. Upon manual examination of
the SD predictions, three primary cases in which failure predictions occur
became apparent: highly noisy CT images, FOVs containing only small portions
of the target structures and pathologies in the target structure (see appendix
fig. 5. While the MD collection may include similar cases, the SDs cover a
less diverse image distribution. As a result, the SD models have larger
variability and perform exceptionally poorly in some rare cases, which may
overshadow the impact of the proposed method. For a further perspective of the
outcomes, a stability ranking analysis is presented in the appendix in fig. 5
[26].
Unfortunately, no corresponding ground truth annotations for the tumor classes
present in the MSD datasets were available in the test dataset, despite the
presence of images with pathologies. Therefore, a clinician carefully
inspected all tumor predictions with a total volume of more than approximately
$1cm^{3}$ to evaluate the effectiveness of our method. The results are
presented in fig. 3 for the MD setting and in table 3 in the appendix. In fig.
3, it can be seen that the proposed region loss successfully mitigates tumor
predictions in anatomical spurious body regions for the MD training setting to
85.7%. Table table 3 in the appendix shows that MultiTalent finds more True
Positives (TP) but also more False Positives (FP) tumors than nnU-Net. Again,
we believe that MultiTalent performs better due to the MD training, but the
Sigmoid activation function allows more FP predictions than the Softmax
function used by default for nnU-Net. Most remaining FP predictions correspond
to anatomical regions within the valid axial region, such as false liver tumor
predictions in the spleen. In the SD setting, the region loss has no
observable influence on reducing the FP predictions, but training with the
unlabeled support set dataset seems to have a positive influence on finding
more true positive tumors in the relevant regions. We will discuss this more
thoroughly in the next section. As it was ambiguous to determine TP and TR for
the Colon MSD task, due to the difficulty of the task itself, we have omitted
this dataset [2].

Figure 3: The proposed Region Loss mitigates most anatomical implausible False
Positive tumor predictions in the MD setting.

##  4 Discussion

The paper presents a novel method aimed at enhancing generalization across
varying FOVs where supervision was initially not feasible using the available
training data. This innovative approach utilizes a BPR model to encode axial
slice scores, enabling the identification of anatomical regions where the
introduced Region Loss function penalizes implausible predictions. Our method
demonstrates superior performance in both MD and SD training settings,
effectively reducing false positive predictions in anatomically unreasonable
regions.
In contrast to post-processing methods, our approach offers greater stability,
as post-processing techniques may inadvertently remove true positive
predictions and rely on an additional error-prone step. Indeed, we encountered
situations where manual correction of the slice scores was required while
applying the BPR model to the images. In general, our results underscore the
robustness of MultiTalent, particularly attributable to its MD training
approach, compared to SD models. While MultiTalent’s use of the Sigmoid
activation function may lead to more false positive predictions, our proposed
method effectively mitigates this issue. The significance of the results in
the SD setting was limited. Including an additional validation set would
enable us to fine-tune our method’s hyperparameters more effectively. However,
it is challenging to find datasets with matching annotated target structures
and an extended FOV. Additionally, because the unlabeled support set came from
a single data source, the models may just learn to ignore predictions for
images within that distribution. As part of our future work, we aim to
integrate a more diverse support set, akin to the MD framework, which is
expected to improve the efficiency of our method. We have evaluated the
efficacy of our method for segmentation in CT images, but it can potentially
be extended to all 3D modalities and tasks such as detection.

## References

  * [1] Antonelli, M., Reinke, A., Bakas, S., Farahani, K., AnnetteKopp-Schneider, Landman, B.A., Litjens, G., Menze, B., Ronneberger, O., Summers, R.M., et al.: The medical segmentation decathlon. arXiv:2106.05735 (2021)
  * [2] Antonelli, M., Reinke, A., Bakas, S., Farahani, K., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., Ronneberger, O., Summers, R.M., et al: The medical segmentation decathlon. Nature Communications (2022)
  * [3] Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., Tarbox, L., Prior, F.: The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository. Journal of Digital Imaging (2013)
  * [4] Daza, L., Pérez, J.C., Arbeláez, P.: Towards robust general medical image segmentation. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2021 (2021)
  * [5] "Gatidis S, a.T.: A whole-body fdg-pet/ct dataset with manually annotated tumor lesions (fdg-pet-ct-lesions). The Cancer Imaging Archive, (2022)
  * [6] Gibson, E., Giganti, F., Hu, Y., Bonmati, E., Bandula, S., Gurusamy, K., Davidson, B., Pereira, S.P., Clarkson, M.J., Barratt, D.C.: Automatic multi-organ segmentation on abdominal ct with dense v-networks. IEEE Transactions on Medical Imaging (2018)
  * [7] He, X., Yang, S., Li, G., Li, H., Chang, H., Yu, Y.: Non-local context encoder: Robust biomedical image segmentation against adversarial attacks. Proceedings of the AAAI Conference on Artificial Intelligence (2019)
  * [8] Heller, N., Sathianathen, N., Kalapara, A., Walczak, E., Moore, K., Kaluzniak, H., Rosenberg, J., Blake, P., Rengel, Z., Oestreich, M., Dean, J., Tradewell, M., Shah, A., Tejpaul, R., Edgerton, Z., Peterson, M., Raza, S., Regmi, S., Papanikolopoulos, N., Weight, C.: The kits19 challenge data: 300 kidney tumor cases with clinical context, ct semantic segmentations, and surgical outcomes. arXiv:1904.00445 (2020)
  * [9] Huang, Z., Wang, H., Deng, Z., Ye, J., Su, Y., Sun, H., He, J., Gu, Y., Gu, L., Zhang, S., Qiao, Y.: Stu-net: Scalable and transferable medical image segmentation models empowered by large-scale supervised pre-training (2023)
  * [10] Isensee, F., Jaeger, P.F., Kohl, S.A.A., Petersen, J., Maier-Hein, K.H.: nnU-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature Methods 18(2)(2), 203–211 (2021)
  * [11] Isensee, F., Ulrich, C., Wald, T., Maier-Hein, K.H.: Extending nnu-net is all you need. In: Bildverarbeitung für die Medizin 2023 (2023)
  * [12] Ji, Z., Guo, D., Wang, P., Yan, K., Lu, L., Xu, M., Zhou, J., Wang, Q., Ge, J., Gao, M., Ye, X., Jin, D.: Continual segment: Towards a single, unified and accessible continual segmentation model of 143 whole-body organs in CT scans. In: IEEE/CVF International Conference on Computer Vision, ICCV 2023 (2023)
  * [13] Kloenne, M., Niehaus, S., Lampe, L., Merola, A., Reinelt, J., Roeder, I., Scherf, N.: Domain-specific cues improve robustness of deep learning-based segmentation of ct volumes. Scientific Reports (2020)
  * [14] Lambert, Z., Petitjean, C., Dubray, B., Ruan, S.: Segthor: Segmentation of thoracic organs at risk in ct images. arXiv:1912.05950 (2019)
  * [15] Landman, B., Xu, Z., Igelsias, J.E., Styner, M., et al.: 2015 miccai multi-atlas labeling beyond the cranial vault workshop and challenge. In: Proc. MICCAI Multi-Atlas Labeling Beyond Cranial Vault—Workshop Challenge (2015)
  * [16] Li, H., Zhou, J., Deng, J., Chen, M.: Automatic structure segmentation for radiotherapy planning challenge (2019), https://structseg2019.grand-challenge.org/ 25/02/2022)
  * [17] Lyu, S., Fan, Y., Ying, Y., Hu, B.G.: Average top-k aggregate loss for supervised learning. IEEE transactions on pattern analysis and machine intelligence 44 (2022)
  * [18] Ma, J., Zhang, Y., Gu, S., Ge, C., Ma, S., Young, A., Zhu, C., Meng, K., Yang, X., Huang, Z., Zhang, F., Liu, W., ant et al., Y.P.: Unleashing the strengths of unlabeled data in pan-cancer abdominal organ quantification: the flare22 challenge. arXiv preprint arXiv:2308.05862 (2023)
  * [19] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. Lecture Notes in Computer Science (2015)
  * [20] Roth, H.R., Lu, L., Farag, A., Shin, H.C., Liu, J., Turkbey, E., Summers, R.M.: Deeporgan: Multi-level deep convolutional networks for automated pancreas segmentation. arXiv:1506.06448 (2015)
  * [21] Roth, H.R., Lu, L., Farag, A., Shin, H.C., Liu, J., Turkbey, E.B., Summers, R.M.: Deeporgan: Multi-level deep convolutional networks for automated pancreas segmentation. In: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015 (2015)
  * [22] Roy, S., Koehler, G., Ulrich, C., Baumgartner, M., Petersen, J., Isensee, F., Jäger, P.F., Maier-Hein, K.H.: Mednext: Transformer-driven scaling of convnets for medical image segmentation. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2023 (2023)
  * [23] Schuhegger, S.: Mic-dkfz/bodypartregression: (v1.0). Zenodo: https://doi.org/10.5281/zenodo.5195341 (2021)
  * [24] Ulrich, C., Isensee, F., Wald, T., Zenk, M., Baumgartner, M., Maier-Hein, K.H.: Multitalent: A multi-dataset approach to medical image segmentation. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2023 (2023)
  * [25] Wasserthal, J., Breit, H.C., Meyer, M., Pradella, M., Hinck, d., Sauter, A.W., Heye, T., Boll, D.T., Cyriac, J., Yang, S., Bach, M., Segeroth, M.: Totalsegmentator: Robust segmentation of 104 anatomic structures in ct images. Radiol Artif Intell. (2023)
  * [26] Wiesenfarth, M., Reinke, A., Landman, B.A., Eisenmann, M., Saiz, L.A., Cardoso, M.J., Maier-Hein, L., Kopp-Schneider, A.: Methods and open-source toolkit for analyzing and visualizing challenge results. Scientific Reports (2021)
  * [27] Yan, K., Lu, L., Summers, R.M.: Unsupervised body part regression via spatially self-ordering convolutional neural networks. In: 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018). IEEE (2018)
  * [28] Zhang, L., Wang, X., Yang, D., Sanford, T., Harmon, S., Turkbey, B., Wood, B.J., Roth, H., Myronenko, A., Xu, D., Xu, Z.: Generalizing deep learning for medical image segmentation to unseen domains via deep stacked transformation. IEEE Transactions on Medical Imaging (2020)

## Appendix

Figure 4: a) An example of unreasonable hepatic liver tumor predictions of
MultiTalent and nnU-Net in the head that our proposed Region Loss completely
mitigates. b) Example spleen predictions for which MultiTalent is more robust
compared to nnU-Net.

Figure 5: Comparison of the Ranking stability of all methods following
Wiesenfarth et al. [26]. The area of each blob is proportional to the relative
frequency of the achieved rank across 1000 bootstrap samples. The median rank
for each algorithm is indicated by a black cross. 95% bootstrap intervals
across bootstrap samples are indicated by black lines. Table 1: List of all
MSD 5-fold cross-validation Dice scores for the MultiTalent networks trained
with different valid boundaries between $\mu$ and $\mu\pm 6\sigma$ to
determine the best boundaries for each class. The networks were trained for
1000 epochs to reduce compute overhead and without a Gaussian boundary
weighting to find the best-fitting boundaries for each class.

class mean slice scores Dice for different offsets to the mean slice scores no
name minimum maximum $0\sigma$ $0.5\sigma$ $1\sigma$ $1.5\sigma$ $2\sigma$
$2.5\sigma$ $3\sigma$ $4\sigma$ $5\sigma$ $6\sigma$ 1 liver $33.09\pm 3.31$
$52.73\pm 1.88$ 95.72 95.80 95.82 95.79 95.89 95.84 95.85 95.84 95.89 95.81 2
liver tumor $39.4\pm 5.46$ $48.49\pm 4.39$ 64.25 62.55 63.55 62.79 63.15 62.82
64.23 61.99 62.13 62.59 3 lung tumor $61.25\pm 4.5$ $64.55\pm 5.15$ 67.23
67.88 65.95 67.54 67.35 65.43 66.14 68.02 68.60 66.68 4 pancreas $34.34\pm
1.96$ $44.1\pm 1.68$ 79.86 79.97 79.97 80.01 80.03 79.85 79.90 79.78 79.81
80.02 5 pancreas cancer $37.67\pm 2.74$ $40.38\pm 2.88$ 55.29 54.51 55.07
55.84 55.62 55.32 55.76 55.11 55.27 55.56 6 hepatic vessel $36.14\pm 4.09$
$51.61\pm 1.74$ 63.76 64.17 64.14 63.92 63.76 64.07 63.89 64.00 64.01 63.89 7
liver cancer $40.74\pm 4.66$ $47.65\pm 3.96$ 70.90 70.95 70.27 70.71 70.92
70.78 70.55 69.75 69.93 70.66 8 spleen $39.87\pm 2.34$ $50.37\pm 1.93$ 96.13
96.07 96.14 96.26 96.29 96.26 96.20 96.23 96.30 96.27 9 colon cancer $24.03\pm
8.83$ $29.49\pm 9.29$ 52.96 52.41 51.91 51.21 53.10 51.31 50.96 52.42 51.04
49.79

Table 2: Overview of the training hyperparameters. For the MD training, the
target spacing, patch size, and number of epochs were fixed as in the original
publication [24]. Additionally, each dataset is sampled with a probability
that is inversely proportional to $\sqrt{n}$, where $n$ is the number of
images per dataset. For all experiments, we selected a batch size of two. For
the SD training, we randomly selected unlabeled images with a 20% probability
and increased the default nnU-Net number of training epochs accordingly.

Framework Setting Batch size Spacing Patch size Sampling prob. epochs $\alpha$
K MultiTalent MD 2 [1.5,1,1] [96,192, 192] default 2000 1 1 nnU-Net MSD Liver
2 [1 ,0.77 , 0.77 ] [128, 128, 128] 0.2 1200 10 1 nnU-Net MSD Lung 2 [1.24,
0.79, 0.79] [80, 192, 160] 0.2 1200 10 1 nnU-Net MSD Pancreas 2 [2.5, 0.80,
0.80] [40, 224, 224] 0.2 1200 10 1 nnU-Net MSD H.Vessel 2 [1.50, 0.80, 0.80]
[64, 192, 192] 0.2 1200 10 1 nnU-Net MSD Spleen 2 [1.60, 0.79, 0.79] [64, 192,
160] 0.2 1200 10 1 nnU-Net MSD Colon 2 [3.0, 0.78, 0.78 ] [56, 192, 192] 0.2
1200 10 1

Table 3: The analysis of all tumor predictions larger than approximately
$1cm^{3}$ on the test set shows the reduction of False Positive (FP)
predictions using the proposed Region Loss. A clinician classified all tumor
predictions to be True Positive (TP) if the prediction certainly indicates a
tumor in the target region (e.g. a tumor within the target region liver).
Predictions within the corresponding target region, where a positive
classification is not unambiguous, for instance, due to fluctuating image
quality, are marked as True Region (TR). Predictions beyond the target region
are always classified as FP.

Method Class TR TP FP MultiTalent Liver Tumor 18 12 23 Lung Tumor 6 13 25
Pancreas Tumor 13 0 19 Hepatic Tumor 14 12 24 Sum 51 37 91 nnU-Net Liver Tumor
17 10 10 Lung Tumor 3 7 8 Pancreas Tumor 9 0 0 Hepatic Tumor 7 6 19 Sum 36 23
37 Method Class TR TP FP MultiTalent + Region Loss (proposed) Liver Tumor 15
12 1 Lung Tumor 3 17 7 Pancreas Tumor 11 0 2 Hepatic Tumor 15 13 3 Sum 44 42
13 nnU-Net + Region Loss (proposed) Liver Tumor 15 12 8 Lung Tumor 3 10 9
Pancreas Tumor 8 0 0 Hepatic Tumor 10 5 19 Sum 36 27 36