

Department of Social Systems and Management
Discussion Paper Series

No.1315

Pooling Experiments for Consecutive Positives
bySong LUO, and Mingchao ZHANG

January 2014

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPAN

1

Pooling Experiments for Consecutive Positives
Song Luo and Mingchao Zhang
Abstract—In this paper, we study the pooling experiments for screening clone maps. Clones are overlapped and placedin a linear order, and hence the clones containing a particular DNA sequence of our interest form a consecutive set (orseveral consecutive subsets) under the linear order. Those clones are called consecutive positives. Pooling experiments areused to identify as many positives as possible by screening as few pools as possible. Stochastic models are introduced to
explain the usefulness of the overlap structure for more efficient pooling experiments even when experimental errors exist. We
describe an efficient positive detecting algorithm, called modified Markov chain pool result decoder (MMCPD), for consecutive
positives. We also introduce an efficient algorithm, called random pooling designer (RPD), for estimating the minimal numberof randomly generated pools required for achieving the complete identifiability of pooling experiments when consecutive
positives and experimental errors exist. To our best knowledge, MMCPD is the first probabilistic group testing algorithm fordetecting consecutive positives when experimental errors exist. Interestingly, RPD shows some happy coincidences between ourtheoretical computation results and previously known simulation results. Some simulation results are also reported in hoping ofdemonstrating that MMCPD and RPD may be promising for real settings.
Index Terms—pooling experiment, group testing, DNA library screening, clone map, consecutive positive, random k-set design,
information-theoretic lower bound.

✦

1

I NTRODUCTIONGROUP testing was proposed by Robert Dorfman[15] during World War II in order to efficiently
test a large number of blood samples for a raredisease. Since then, various applications of group
testing have been found in many fields such as multiple communication, coding theory, information security, sparse signal recovery and others. Particularly,
biology-motivated group testing has been developed
into one of the most important tools in the studyof gene functions. For example, in Human Genome
Project, well-designed group testing schemes havebeen proved to be useful in both saving materials andaccelerating the process of reconstructing high-quality
DNA libraries. These high-quality libraries have beenfrequently and repeatedly used for extensive studies.
A DNA library is a collection of cloned DNA segments taken from a specific organism. Those DNA
segments are called clones. Determining whether aclone contains a particular DNA sequence of interest
can be accomplished by screening it with a probe.
The clone is called a positive for the probe if it
contains the particular DNA sequence, and a negative
otherwise. Due to the large size of a DNA library andthe cost of materials, instead of screening each clone
individually, combinations of the clones are screened
to identify and isolate the positives. The combinationsof the clones are called pools. Each pool is screened• S. Luo is with the Graduate school of Systems and Information
Engineering, University of Tsukuba, Ibaraki, Japan, 305-8577.
E-mail: luo@sk.tsukuba.ac.jp• M. Zhang is with the Graduate school of Systems and Information
Engineering, University of Tsukuba, Ibaraki, Japan, 305-8577.
E-mail: mzhang@sk.tsukuba.ac.jp

with the probe to learn whether any of the clones in
the pool contain the DNA sequence. Screening poolsin this way is called a pooling experiment. Poolingexperiments are used to identify as many positives
as possible by screening as few pools as possible. The
efficiency of pooling experiments has been studied byBarillot et al. [2], Berger et al. [3], Bruno et al. [5] andSham et al. [38]. Here we refer to the books [12] and[13] by Du and Hwang for an overview.
1.1

Classification of Group Testing

Any group testing consists of a pooling procedure anda positive detecting procedure. A pooling procedure isa procedure of constructing a collection of pools called
pooling design, determining which clones are put
into which pools. A positive detecting procedure is aprocedure of determining which clones are positives
from the outcomes of group tests.
Group testing can be classified as either adaptive or
nonadaptive, based on how a pooling design is constructed in the pooling procedure. In adaptive group
testing, a pooling design is constructed in multiplestage. In each stage, the outcomes of the group tests
from previous stages are learned to construct thepools in the next stage. In nonadaptive group testing,a pooling design is constructed in one-stage, before
any outcome of the group test is known. In practice,
it is preferable to screen pools in parallel or with few
stages for minimizing time consumption.
Group testing can also be classified as either combinatorial or probabilistic.To implement a combinatorialgroup testing, we have to construct a pooling design
with desirable combinatorial properties such that all
positives can be distinguished from negatives under2

the assumption on the maximum number of positives
and that of experimental errors. Related studies can befound in Du and Hwang [12], [13], [14], Dyachkov etal. [16] , Macula [29], and Ngo and Du [33]. To implement a probabilistic group testing, we need not only
stochastic models for positives and for pooling results
but also an efficient positive detecting algorithm toinfer positives from erroneous pooling results basedon a stochastic inference model. Probabilistic group
testing is developed by Bruno et al. [5], Knill et al. [27],Mezard and Toninelli [30] and Uehara and Jimbo [37].
Bruno et al. [5] and Knill et al. [27] proposed a positive
detecting algorithm called MCPD by using Markov
chain Monte Carlo simulation method. Uehara andJimbo [37] proposed another efficient algorithm called
BNPD by using Bayesian network inference. As far as
we know, MCPD and BNPD are the only known efficient positive detecting algorithms when experimental
errors exist.
Probabilistic group testing has three noticeable merits in practice. First, from the perspective of detecting
procedure, probabilistic group testing is expected toperform stably even when a relatively larger numberof positives or/and experimental errors occur than
expected. Second, probabilistic group testing reduces
information loss by interpreting a measurement of atest into a multilevel state. Third, from the perspective of pooling procedure, implementing probabilistic
group testing requires few restriction on the combinatorial structure of pooling designs. Particularly,
random pooling designs are allowed. In fact, random
pooling designs are often preferred in a real setting.
This is because (i) a well-designed random poolingdesign can have a satisfying efficiency with desirable
error-tolerant ability, (ii) it is unrealistic to expect tobe able to find an appropriate combinatorial poolingdesign for every new pooling experiment, and (iii)
random pooling designs facilitate robot automationand hence are very easy to construct.
1.2

Linear DNA library and Consecutive Positives

Motivated from applications to DNA library screening, Balding and Torney [4] considered the problemof pooling experiments for screening unique-sequence
on a 1530-clone map of Aspergillus nidulans. The
clone map has the properties that the clones are,with possibly a few exceptions, linearly ordered andno more than two of them cover any point on the
genome. The goal of screening clone maps is to identify where a particular DNA segment occurs on aclone map.
In this problem, the clones are overlapped, andhence it may happen that one segment of interest
occurs in a relatively large number of clones, but
typically the number is predictable as it is related tothe clone coverage. Colbourn [9] introduced the dconsecutive property to capture the overlap structurewith efficient combinatorial group testing algorithms.
Following his work, related studies can be found,for example, in Müller and Jimbo [31] and[32], Juanand Chang [24] and Ge et al. [18]. Bruno et al. [6]used heterogeneous priors and discussed optimization issues of nonadaptive random pooling designs on
the assumption that an effective detecting procedure
exists and experimental errors do not exist.
However, to our best knowledge, when experimental errors exist, probabilistic group testing for consecutive positives has not yet been discussed. In this paper,
we study this problem to fill the gap.

2

S TOCHASTIC M ODELSThis section introduces a prior probability distribution for consecutive positives with overlap structure,and extends Knill et al [27]’s method for consecutive
positives.
2.1

Prior Knowledge of Consecutive Positives

The following notation will be consistently usedthroughout. Let C = {c1 , c2 , ..., cn } be a set of clones.
C is called a DNA library. Each clone ci has an
associated state σi ∈ {0, 1}. ci is a positive if σi = 1,otherwise a negative. Let P = {ci : σi = 1} andXP = (σ1 , . . . , σn ) be the set of positives and its
vector form corresponding to P , respectively. For convenience, P and XP are used interchangeably, referredto the positive set. Denote by xi the random variableof σi such that
0, if σi = 0,xi =1, if σi = 1,

and by X = (x1 , . . . , xn ) the random vector of theassociated states of C.
When a DNA library is constructed by uniformly
and randomly cloning of DNA segments, the expectedn���number of positives E
xi = d is used as the prior
i=1

knowledge, for some positive number d. This priorknowledge appeared in Knill [26], Knill et al. [27] andUehara and Jimbo [37]. The overlap information cananalogously be collected and expressed in a similar
form. The basic idea is to extend the d-consecutive
positives property proposed by Colbourn [9] into a
probabilistic version with broader sense.
C is said to be linear if it is associated with an linear
order ci ≺ ci+1 , for 1 ≤ i < n. The positive set of linear
C is said to have the multi-d-consecutive property if
any subset of P that forms a consecutive set (under
the ordering ≺) contains at most d positives. A subset
of positives P � is said to be a maximum consecutive
subset, if P � is a consecutive set but P � ∪ {c} is nota consecutive set for any c ∈ P \ P � . Notice that
P corresponds to a unique partition with maximum

3

consecutive subsets. The positive set with the multid-consecutive property is allowed to have more than
one maximum consecutive subsets, each of whichcontains at most d consecutive positives.
Without loss of generality and for the sake of simplicity, we analyze the overlap structure of the clone
map of Aspergillus. When the clone map is screened,
if the positive set is believed to approximately havethe multi-2-consecutive positive property, then thefollowing prior knowledge of the positives can becollected:
Info 1: the portion of positives among all clones isrelatively small;
Info 2: although the positives are sparse, some of thepositives tend to be consecutive;
Info 3: although some of the positives tend to beconsecutive, maximum consecutive subsets
tend to not get close to each other;
Info 4: no more than two of the positives tend toform a maximum consecutive subset.
Therefore, based on Info 1 through Info 4, we can
express the prior knowledge of consecutive positives
as follows:
























E
E
E
E

n
��

i=1

i=1
� n−2
�
i=1
� n−2
�

�
�

(I)

xi xi+2 = d3
�σj : ∃t ∈ N and ∃ t distinctive sets
C

Since Shannon’s theorem [39] established the uniqueness of entropy as an information measure of uncertainty, the principle of maximum entropy has been
wildly used to derive prior probability distributions.
Intuitively speaking, more information means less uncertainty. Any probability distribution satisfying theconstraints that has less uncertainty will contain more
information, and thus implies something stronger
than what the prior knowledge means. The principle of maximum entropy, as a method of statisticalinference, is due to Jaynes [21], [22] and [23].
Based on the principle of maximum entropy, toderive the prior probability distribution P(X) which
incorporates (I) but is free from any other knowledge,
is to solve the following problem,
−









subject to

C1 , . . . , Ct ∈ 2 \ {∅} such that C ⊆

t
�

i=1

�

Ci .

Each polynomial of D is called an overlap polynomial of C. We can collect and explicitly express priorknowledge of consecutive positives by using overlap
polynomials in accordance with overlap information.
If s overlap polynomials Dj� have been chosen from D,
P is said to have positive pattern d if (dd)j = Dj� (XP ),for j = 1, . . . , s. In addition, we denote by |dd| thenumber of positives in positive pattern d .

�P(X) log P(X)X∈{0,1}n









xi xi+1 xi+2 = d4

Generally, we denote by D a family of polynomialsof σ1 , . . . , σn , wherei=1 cj ∈Ci

The Principle of Maximum Entropy

xi xi+1 = d2

Assumption 1. Appropriate values of d1 through d4 havebeen given as parameters.

D=

2.2.1

P(X)

for some positive numbers d1 , d2 , d3 and d4 . However,
how to assign proper values to d1 through d4 is notwithin our present concern. Formally, we introduce
the following assumption.

��t�This part seeks a prior probability distribution that, insome sense of optimality, approximately incorporates
the prior knowledge of consecutive positives (I).

maximize

xi = d1

� n−1
�

i=1

�

2.2 Prior Probability Distribution of Consecutive
Positives��
EP D1 (X) = d1
��
EP D2 (X) = d2
��
EP D3 (X) = d3
��
EP D4 (X) = d4
�X∈{0,1}n P(X) = 1.

(II)

As is well-known, Lagrange multiplier method leads
to the solutionPθ (X) =

1
exp
Z(θ)

�

−

4
�
j=1

�
θj Dj (X) ,(III)

for some constant vector θ = (θ1 , θ2 , θ3 , θ4 ) and constant Z(θ). In literature, (III) is called a Gibbs measure.
The normalization constantZ(θ) =

�X∈{0,1}n

exp

�

−

4
�

θj Dj (X)j=1

�is called the partition function. It connects θ with theconstants di by simultaneous equations∂ − log Z(θ)= di ,∂θi(IV)for i = 1, . . . , 4. Before discussing how to derive θfrom (IV), we first point out three useful properties of(III).
n���Property 1 (Consistency Property). If E
xi = di=1(d > 0) is the only prior knowledge of positives, then the

4prior probability distribution determined by the principleof maximum entropy is
n
n
� ��
�n− �
x �
xi
d i=1 id
i=1
P(X) =
1−.
nn

The proof can be found in Jaynes [23]. From thisproperty, we see that, provided proper prior knowledge of positives, the principle of maximum entropy
can be used to obtain the prior probability used in
Bruno et al [5], Knill [27], Knill et al [26] and Uehara[37], where the DNA library is considered to havebeen constructed by uniformly and randomly cloningof DNA segments. Hence, maximum-entropy distributions may be reasonable extensions for describing
the overlap structure of clone maps.
Denote by Nj the set of neighbors of j. Define N1 �
{2, 3}, N2 � {1, 3, 4}, Nn−1 � {n − 3, n − 2, n}, Nn �{n − 2, n − 1} and Ni � {i − 2, i − 1, i + 1, i + 2}, for3 ≤ i ≤ n−2. For any configuration of xi by assigning
xi = σi for each i ∈ [n] according to Pθ (X), (III) has
the following Markov-type property.
Property 2 (Local Markov Property). For any j ∈ [n],Pθ (xj |σi , i ∈ [n] \ {j}) = Pθ (xj |σi , i ∈ Nj ).
We show a sketch of the proof. A more general one
can be found in Kindermann and Snell [25] and Pearl[34]. To see this, we rewrite the Gibbs measure into
the product form:Pθ (X) =

n
n−1
�
1 �
exp{−θ1 xi }exp{−θ2 xi xi+1 }Z(θ) i=1
i=1
n−2
�
i=1

exp{−θ3 xi xi+2 }

n−2
�
i=1

exp{−θ4 xi xi+1 xi+2 }.

The proof is obvious. Notice that, for any fixed
positive pattern d , Pθ is the heterogenous prior withthe greatest uncertainty.
2.2.2

The constants θj are called Lagrange multipliers. As
the same values of dj ’s may come from different information sources, different estimation methods may
be needed to derive θ from (IV). Here we discuss twocases. In the first case, the values of dj ’s represent thedegrees of belief. The belief is needless to be relevantwith the outcomes of any random experiment. Whilein the second case, the values are obtained from theobservation of the data Y = {Y1 , . . . , YN }.
In the former case, we definef (θ) = −4�j=1

dj θj − log Z(θ).

From (IV), we see that θ can be any stationary pointof f (θ). Therefore, to obtain θ is to solve a gradientsquare-minimization problem.
minimize �∇f (θ)�2 .
θ∈R4

log

N
�

i=1

P(Yi |θ) =

N �
�
i=1

−

Property 3 (Heterogenous Property). Given X1 , X2 ∈{0, 1}n , if they have the same positive pattern, thenPθ (X1 ) = Pθ (X2 ).

4
�
j=1

�

Nj (Yi )θj − log Z(θ) .
(∗)

Let lY (θ) = N1 log P({Y1 , . . . , YN }|θ) and dj =N�1
Dj (Yi ), for j = 1, . . . , 4. To obtain a maximumN
i=1

likelihood estimate of θ is to maximize (∗). Equivalently, we solve

σj =0Thus, after canceling out the normalization constant,terms of Pθ (σ1 , . . . , σn ) that do not contain σj cancel
from both the numerator and denominator of thecondition probability and therefore this probability
depends only on the random value of xj and thoseof its neighbors. This property serves a probabilistic
interpretation of multi-2-consecutive property, sayingthat the state of a clone is only influenced by the statesof its previous two and also next two neighbors, if theclone has such neighbors.

(V)

However, in the latter case, we may fit a model
chosen from {Pθ : θ ∈ R4 } to the given data.
Assuming that Y1 , . . . , YN are i.i.d (independent andidentically distributed) random samples drawn fromPθ with unknown parameter θ, we estimate θ by applying maximum likelihood estimation method. The
log-likelihood is defined by

For σj = 0, 1, the condition probability on the left
side is,Pθ (xj = σj |σi , i ∈ [n] \ {j})Pθ (σ1 , . . . , σj−1 , xj = σj , σj+1 , . . . , σn )= 1.
�Pθ (σ1 , . . . , σj−1 , xj = σj , σj+1 , . . . , σn )Estimation of Lagrange Multipliermaximize lY (θ) = −
θ∈R4

4
�
j=1

dj θj − log Z(θ).

(VI)

Both problems are very difficult to solve, because
the partition function Z(θ) involves exponentiallymany computations and usually cannot be known
thus. Without calculating the partition function, Geyer
and Thompson [20] developed a method to numerically solve (VI) by using importance sampling andMonte Carlo simulation. Descombes et al [11] further demonstrated an efficient conjugate gradient algorithm. Motived by their work, we employ their
methods to solve (V). Here we sketch that (V) is alsosolvable (see [20] and [11] for detailed discussions).

5

The key idea is to estimate, for any fixed θ,�∇f (θ)�2 and its gradient by using importance sampling. By employing
� �
�
4
1
Pψ (X) =exp −ψj Dj (X) ,Z(ψ)
j=1
we reformulate f (θ) as follows:f (θ) = −It follows that

Z(θ) =

�

4
�
j=1

exp

X

· exp

�

�

dj θj − log�−= Eψ exp

−

4
�
j=1

4
�

Z(θ)
− log Z(ψ).
Z(ψ)(θi − ψi )Dj (X))ψj Dj (X)

i=j

�

−

4
�
j=1

�

maximize f (θ) = −

�

�

θ∈R4

�

(θi − ψi )Dj (X) Z(ψ) ,where Eψ refers to the expectation with respect to Pψ .
Then, we obtain
�
� �
��
4
Z(θ)
= Eψ exp −(θj − ψj )Dj (X) .
Z(ψ)
j=1

(2.1)

The significance of (2.1) lies in that although Z(θ)Z(θ)
is unknown, Z(ψ)
can be estimated by a sampling of
the known probability distribution Pψ . Furthermore,
1 ∂Z(θ)we can rewrite Z(ψ)∂θj into�� ���
4Eψ − Dj (X) exp −(θi − ψi )Dj (X) ,

(2.2)

j=1

2

1 ∂ Z(θ)
and similarly Z(ψ)∂θj ∂θk into�Eψ Dj (X)Dk (X) exp

�

−

4
�
j=1

(θi − ψi )Dj (X)(2.1), (2.2), (2.3) and even higher-order partial deriva1
tives of Z(θ) up to constant Z(ψ). This allows us
to learn local variation of �∇f (θ)�2 at any fixed θ.
Those estimations will be helpful to numerically sovle(V). Particularly, the local markov property of Pψ
makes the sampling procedure efficient. We refer to
Descombes et al [11] for other algorithmic details.
In fact, (VI) is stronger than (V) in the sense thatany optimal of (VI) is also an optimal of (V), but not
the reverse. Therefore, we may also formulate the firstcase in the stronger sense,

��

.

(2.3)

Next, with (2.1), (2.2) and (2.3) we can obtain,∂�∇f (θ)�2
∂θk
�� �
�2 �
4 �
�
Z(ψ) 1 ∂Z(θ)
Z(ψ)
1
=2
dj +
−
Z(θ) Z(ψ) ∂θj
Z(θ)
Z(ψ)
j=1
��
�
�
��
2
∂Z(θ)
1 ∂Z(θ)Z(ψ)
1 ∂ Z(θ)+.
∂θk
Z(ψ) ∂θj
Z(θ) Z(ψ) ∂θj ∂θk(2.4)
Using the Gibbs sampler (Geman and Geman [19]),
for any fixed θ, by Pψ we can theoretically estimate

4
�
j=1

dj θj − log Z(θ).

(VII)

From this formulation, we can see that both cases
may be numerically solvable, but the first case ismore general than the second one. By estimating a
numerical solution of (V) within a given tolerance, themodel can be approximately determined, and it canbe served as the desired prior probability distributionthat approximately incorporates our prior knowledge
derived from the available information.
Assumption 2. Given the values of d1 through d4 under
Assumption 1, θ� is an approximate optimal of (VII). Pθ�(X)is the desired prior probability distribution.
2.3

Stochastic Model of Pooling Results

The following notation concerning pooling experiments will be consistently used. Let A = {p1 , . . . , pm }be a pooling design. Each pool pi is a subset of C
corresponding to the clones in the pool. The poolingdesign A constructed by determining each of n clones
is put into which of m pools, can be represented byan m × n binary matrix A = (aij ). Each entry aij isdefined as follows:�
1, if cj ∈ pi ,aij =0, if cj �∈ pi .
Then, A is called the incidence matrix of poolingdesign A and m the size of A. For convenience, A and
A are interchangeably referred to a pooling design.
The pooling result of screening pool pi is denotedby r(pi ). An actual observation of r(pi ) is often given
in multilevel measurement, taking any value from a
set of test outcomes V , such as ”negative”, ”weakpositive”, ”medium positive” or ”strong positive”.
Unfortunately, the pooling results are typically corrupted due to the inclusion of additional clones in apool or to the failure of screening test. Therefore, weneed infer P from erroneous pooling results.
To begin with, we introduce the stochastic model ofpooling results given in Knill et al. [27]. The following
assumptions will be used.

6

Assumption 3 (Knill et al. [27]). The distribution ofr(pi ) depends only on the number of positives, |pi ∩ P |,in pi .
Assumption 4 (Knill et al. [27]). Since the PCR was used
for screening, we assume that Pr(r(pi )||pi ∩ P |) depends
only on |pi ∩ P | = 0 or |pi ∩ P | ≥ 1.
For any integer b, we define�
0 if b = 0,b=
1 if b ≥ 1.

Notice that, for each i, the number of positives inpool pi can be represented by the ith coordinate of
AXP , that is, |pi ∩ P | = (AXP )i . Let Ri be a randomvariable of pooling result r(pi ),
0 if r(pi ) is negative,
 1 if r(p ) is weak positive,iRi =2 if r(pi ) is medium positive,
3 if r(pi ) is strong positive.
Rewrite Pr(Ri = ri |(AXP )i ) into f (ri , (AXP )i ) for
short, and denote by PrA (R = r|XP ) the likelihoodfor the pooling results r(∈ {0, 1, 2, 3}m ) obtained frompooling design A with m pools. By using Assumptions 3 and 4 and the chain rule, it can be expressed
as a product:

PrA (R = r|XP ) =m�Assumption 6. Given any A ∈ A, PrA (X) = Pθ�(X),for all X ∈ {0, 1}n .
From Assumption 2, Assumption 5, Assumption 6
and (VIII) , it follows that the posterior probability iscomputable up to a multiplicative constant, and canbe written as,
PrA (XP |r) ∝ Pθ�(XP )m�f (ri , (AXP )i ).

i=1

If we are interest in PrA (c ∈ P |r) for each clone c,this probability can be written as the marginal

PrA (xj = 1|r) ∝

�

P ⊆C:
cj ∈P

Pθ�(XP )m�

i=1

f (ri , (AXP )i ).

(IX)

Notice that, if we substitute Binomial distribution
for Pθ�, the inference model turns out to be the oneused in Knill et al. [27].

3f (ri , (AXP )i ).

(VIII)

i=1

Let f = {f (j, i) : j = 0, . . . , 3, and i = 0, 1}. In realapplications, f can be appropriately estimated (seeKnill. et al [27]). In this application, we assume thatf have been estimated and provided as parameters.
Assumption 5. The values of f (j, i) are known a priori,
but f (j, i) �= 0, for j = 0, . . . , 3, and i = 0, 1.
Particularly, notice that f (1, 0) + f (2, 0) + f (3, 0)
is the likelihood of a false positive, whereas f (0, 1)
is that of a false negative. Moreover, due to Assumption 3 and Assumption 4, the model does not
require any restriction on the structure of positives,
and thus is compatible with the overlap structure ofconsecutive positives.
2.4

to use. Prior knowledge of positives will affect theconstruction of pooling design A. However, the reverse does not reasonably hold, because knowing
the construction of pooling design A will not add
any information to the prior knowledge of positives.
Therefore, we introduce the following assumption.

Bayes Inference Model

To infer the consecutive positives, Bayes inference
model is used to decode the erroneous pooling results.
Denoting by PrA (XP |r) the posterior probability that
XP is the positive set given the pooling results r under
the pooling design A, by Bayes’ rule we have
PrA (XP |r) ∝ PrA (XP )PrA (r|XP ),where PrA (XP ) denotes the probability that XP is
the positive set given pooling design A is chosen

D ETECTING A LGORITHM FOR C ONSECU TIVE P OSITIVESThis section presents a positive detecting algorithm
for consecutive positives, called modified Markov
chain Monte Carlo pool result decoder (MMCPD).
Our primary goal is to compute PrA (xj = 1|r),for each j, and to choose the clones with highest
posterior probabilities of being a positive for confirmatory individual tests. Those clones are calledcandidate positives. However, as it stands in (IX),
every exact computation will involve exponential
work in n, making our primary goal impracticalas n is very large. Instead, we estimate the posterior probabilities by drawing samples approximately
according to PrA (X|r). This idea works as longas a sampling method exists and is able to efficiently produce sufficiently many samples according
to PrA (X|r). To sample from PrA (X|r), we can use the
Gibbs sampler, since the full conditional distributions
PrA (xi |σ1 , . . . , σi−1 , σi+1 , . . . , σn , r) are easy to obtainfor each i, due to the local markov property of Pθ�(X)and to the product expression of Pr(r|X).
In this way, Knill et al. [27]’s method can be extended for detecting consecutive positives. Mimicking their approach, we construct a Markov chain
X0 , X1 , . . . on the family of all configurations of C withPrA (X|r) as its stationary distribution. Denote by XS
the vector form of the positive set S. Let XS be theconfiguration of C at the end of step t − 1, that is,7state Xt−1 . Step t of the chain updates XS by making
a random decision for each of the n clones on whether
to add it to or remove it from S. At the end of step t,that is, after all clones have been processed, state Xt
is the final XS . In the Gibbs sampler, for clone ck , theprobability of changing from XS to XS�{ck } equals to
1,
PrA (XS |r)1+
PrA (XS�{ck } |r)(X)

where � is the symmetric difference, that is, A�B �(A \ B) ∪ (B \ A). Notice that with this probability ofchanging from XS to XS�{ck } , Xt and hence XS tendin distribution to PrA (X|r). From Bayes’ Theorem,Pθ�(X) and (VIII) we obtain,
PrA (XS |r)
PrA (XS�{ck } |r)=�f (ri , (AXS )i )

i∈[m]:Pθ�(XS )ck ∈pi.
�Pθ�(XS�{ck } )f (ri , (AXS�{ck } )i )
i∈[m]:
ck ∈pi

Notice that, for i such that ck ∈ pi , we haveaccording to PrA (X|r). By taking sufficiently many
samples from the Markov chain, the proportion of thesamples with σi = 1 can be thus a Bayesian estimateof PrA (xi = 1|r).
As an extension of MCPD, MMCPD has two attractive merits. First, only local neighborhood relations
and state values are needed to update the prior ratio
and likelihood ratio, making the computation efficient. Second, its implementation requires no explicit
restriction on the structure of pooling designs. This
provides us a considerable freedom to choose andoptimize the pooling designs we want to use, rather
than the ones we have to use. Particularly, MMCPD
is able to decode the pooling results obtained byscreening random pooling designs.

4

R ANDOM k- SET D ESIGNThis section discusses random pooling designs in thepresence of consecutive positives and experimental
errors.
4.1

Motivation and Related work

Among all, we are particularly interested in random
k-set designs. A is said to be a random k-set design if�each column of A is a k-subset of [m] that is generated(AXS )i + 1, if ck �∈ S,(AXS�{ck } )i =independently and uniformly at random. k is called(AXS )i − 1, if ck ∈ S.
the replication number, deciding how many poolsPθ�(XS )The prior ratio P �(XS�{c } ) can also be efficiently each clone will be put into. Since there may exist
θkother positive detecting algorithms for consecutive
updated and we will see that it only depends on thepositives, a practical lower bound of the optimal
associated states of ck ’s neighbors. For convenience,value of m, independent of any positive detectingS
S
let XS = (σ1S , . . . , σkS , . . . , σnS ) and σ−1= σ0S = σn+1=
algorithm, will be attractive. This motivated us to seekS
σn+2
= 0, for any S. If ck �∈ S, thena nontrivial information-theoretic lower bound of m.
� S
Information-theoretic bounds have been extensivelyPθ�(XS )SSS= exp θ�1 + θ�2 σk−1
+ θ�2 σk+1
+ θ�3 σk−2
studied in many fields. For example, related studiesPθ�(XS�{ck } )�can be found in Atia and Saligrama [1], Chan et al
S
S
S
S
S
S
S
+θ�3 σk+2+ θ�4 σk−2
σk−1
+ θ�4 σk−1
σk+1
+ θ�4 σk+1
σk+2.
[7] and [8], Wang et al. [40] and Wainwright [41].
Our computation method shares some commonplaces
If ck ∈ S, theneither in purpose or in formulation with prior work�Pθ�(XS )S
SSsuch as Bruno et al. [6], Knill et al. [26], Wang et al. [40]= exp − θ�1 − θ�2 σk−1
− θ�2 σk+1
− θ�3 σk−2Pθ�(XS�{ck } )� and Wainwright [41]. However, the algorithmic conS
S
S
S
S
S
−θ�3 σk+2− θ�4 σk−2 σk−1
− θ�4 σk−1
σk+1
− θ�4 σk+1
σk+2
. sideration of consecutive positives and experimental
errors with random k-set designs makes our method
From Assumption 4 and that Pθ� is strictly posdifferent from any of them.
itive (that is, Pθ�(X) > 0 for any X), we can seethat the Markov chain is aperiodic, since it has apositive probability of remaining in the same state, 4.2 Problem Formulationand that the Markov chain is also irreducible, since Notice that the Gibbs sampler (Geman and Geman
it is possible to go from any state to any other state. [19]) allows us to estimate the distribution of the posThis assures uniqueness of and convergence to the itive patterns of P �. Hence, to choose proper values ofθstationary distribution PrA (X|r) as well as the ergodic m and k for P � and f , we begin by fixing a positiveθproperty. For discussions of the Gibbs sampler, we re- pattern d and formulate the subproblem in terms of dfer to Roberts and Smith [35] and Tierney [36]. Hence, and f .
n���starting with any state and running Gibbs sampler
Let Ω(dd) ⊆ (σ1 , . . . , σn ) ∈ {0, 1}n :σi = |dd|
algorithm after a suitable warmup period, samples
i=1
obtained from a realization of the Markov chain can be the nonempty subset consisting of all the vectors
approximately be regarded as the desired ones drawn with positive pattern d . If XP belongs to Ω(dd), due to8

the heterogenous property of Pθ�, we think of XP randomly and uniformly distributed over Ω(dd). Denoteby Am,k the set of all k-set designs with size m. Similar
to Assumption 6, only knowing the k-set design does
not change the distribution of XP over Ω(dd). We
formally states this by introducing Assumption 7.
Without causing confusion, when d is fixed andknown, we denote by Pr(X) the probability that X ischosen from Ω(dd) as the positive set and by PrA (X)the probability that X is the positive set given k-set
design A.
Assumption 7. When d is given and the values of m andk have been decided, given any k-set design A ∈ Am,k ,PrA (X) = Pr(X), for all X ∈ Ω(dd).

Given any instance of random k-set design A, apositive detecting algorithm φ with respect to A is amapping from the m-vector observation r to an estimated vector of positives, say of the form XP� = φA (r).
Accordingly, based on the k-set design A, the average
probability of decoding error of any positive detecting
algorithm is defined as�
1perr (A) =Pr[φA (r) �= X|X].
|Ω(dd)|d)X∈Ω(dWe apply Fano’s lemma [10] to lower bound theaverage probability of decoding error. Notice thatfA
X −→ AX −→ r forms a Markov chain, we can arriveat the form used in Wang et al. [40]The average probability of decoding error over Am,kcan be lower bounded as follows:
EA
EA [perr (A)] ≥1 −i=1

�HA (ri ) − EA [HA (r|X)] + 1
log |Ω(dd)|

4.3

EA [HA (rh )] ≤ −where
fd ,m,k (j) =

1
�
i=0

3
�

fd ,m,k (j) log fd ,m,k (j),

j=0

i
λd1−i
,m,k (1 − λd ,m,k ) f (j, i),d||dLemma 2. Pm,k (t) can be iteratively computed, for t =k, . . . , |dd|k.
d||dPm,k (t) =

k
�

w=0

where

d|−1
|d

µm,k,t (w)Pm,k (t − k + w),µm,k,t (w) = δ(t − k + w, k)�t−k+w��m−(t−k+w)�
w

and in particular�m� k−w,k

Using Lemma 2, we can obtain a closed expressionof EA [HA (r|X)].
Lemma 3.
EA [HA (r|X)] =

.

(XI)

A lower bound of EA [perr (A)]Provided n, d and f are fixed and known, we begin
by stating a set of lemmas and a necessary condition
on m and k for EA [perr (A)] = 0, on the conditions thatthe positive set X with positive pattern d is randomly
and uniformly distributed over Ω(dd) and A is random
k-set design, independently and uniformly from Am,k .
Their proofs will be given in appendices.
Lemma 1. For any h ∈ [m],k |dd|) .
m
We introduce some notations for an exact computation of EA [HA (rh |X)]. Let |AX| = |{h ∈ [m] :(AX)h = 1}|. Given an X with positive pattern d ,d||ddenote by PX,m,k (t) the probability that a k-set design
A, uniformly chosen at random from Am,k , satisfies|AX| = t. Recall that the columns of random kset design are independent and uniformly chosen at
d|
d|
|d
|d
random. This suggests that PX,m,k (t) = PX � ,m,k (t) ford||dany X and X � with d. Thus, we can write Pm,k (t) for
short. For any integers a and b, we define an indicate
function of a and b as follows:�
0 if a < b,
δ(a, b) =1 if a ≥ b.
λd ,m,k = (1 −1Pm,k(k) = 1.

HA (r) − HA (r|X) + 1
perr (A) ≥ 1 −.
log |Ω(dd)|

� m
�

and

d|k
|d
�
t=k

where

�
�
d|
|d
Pm,k (t) tHf (1) + (m − t)Hf (0) ,Hf (i) = −for i = 0, 1.

3�f (j, i) log f (j, i),j=0Respectively substituting the results of Lemma 1
and Lemma 3 for EA [HA (r)] and EA [HA (r|X)] in (XI),
we can easily compute an information-theoretic lower
bound of the average probability of decoding error,that is,EA [perr (A)] ≥ 1 −where
I(dd, f, m, k) = −m−3
�I(dd, f, m, k) + 1,
log |Ω(dd)|(XII)

fd ,m,k (j) log fd ,m,k (j)

j=0

d|k
|d
�
t=k

�
�
d|
|d
Pm,k (t) tHf (1) + (m − t)Hf (0) .

9f (0, 0) = 0.871
f (1, 0) = 0.016
f (2, 0) = 0.035
f (3, 0) = 0.078

f (0, 1) = 0.05
f (1, 1) = 0.11
f (2, 1) = 0.27
f (3, 1) = 0.57

f (0, 0) = 0.856
f (1, 0) = 0.126
f (2, 0) = 0.016
f (3, 0) = 0.002

f (0, 1) = 0.02
f (1, 1) = 0.155
f (2, 1) = 0.288
f (3, 1) = 0.537

By fixing the value of m (m = 47, 60, 131), Fig. 1
shows the variation of ILB(|Ω(dd)|, f, m, k) as the values of f and k vary.
Similar to the LDPC (Low-Density Parity-Check)codes (see Gallager [17] and Mackay and Neal [28]),
random k-set designs demonstrates a degree of errortolerant ability, which is closely related with the valueof k. The U-shaped curve of k is in accordance with
our expectation, since either too small or too large
a value of k will weaken the positive detectabilityof random k-set designs. As Uehara and Jimbo [37]mentioned, we also observe that there is a big gap

0.5
0.2

E(Perr)

0.0
−0.1

0.0

−0.2

−0.1

4

6

8

10

12

14

2

4

6

8

10

Replication Number k

Replication Number k(a) f1(b) f2

12

14

Fig. 1: Comparison of ILB(|Ω(dd)|, f, m, k) subject to
f1 and f2
between the values of ILB(|Ω(dd)|, f, m, k) when k = 3
and k = 4.
Comparing Fig. 1 (a) with Fig. 1 (b), the influence off is also noticeable. Uehara and Jimbo [37] proposed
f2 for repairing the unnatural monotonicity of f1 inwhich f (1, 0) < f (2, 0) < f (3, 0). They also simulated the positive detectability of BNPD and MCPD
subject to f2 and showed that the seemingly slight
modification remarkably improves the performance of
MCPD and BNPD. Interestingly, as is shown in Fig. 1,
the variation of ILB(|Ω(dd)|, f, m, k) implies that f1 ismore difficult to cope with, which also coincides with
their simulation results.
4.4

Random Pooling Designer

Given a positive pattern d it is often over complicated to obtain the exact |Ω(dd)|. Instead, we estimate
a lower bound of |Ω(dd)|. We consider the positive
pattern of the form d = (d1 , d2 , 0, 0) for some positive
integer d1 and nonnegative integer d2 . We claim that
−d2
�
� d1�
2
|Ω(dd)| = d1d−d(n − d1 + d2 − 3i + 4) is a lower2
i=1

TABLE 2: f2 : Uehara and Jimbo [37]

0.1

0.2

E(Perr)

0.3

0.3

m=131

−0.2

log |Ω(dd)| ≤ I(dd, f, m, k) + 1.

f Uehara, et. al.

m=47
m=60
m=131

0.4

0.4

m=60

2

TABLE 1: f1 : Knill et al. [27]n=1298 d=4

f Knill, et. al.

m=47

Theorem 4 (Necessity Theorem). EA [perr (A)] vanishes
to 0 only if m and k satisfy the condition

In each subproblem where d is fixed and known,
the Necessity Theorem implies that roughly m =O(log |Ω(dd)|) pools are least required for complete
identifiability. This casts a light on the usefulness ofknowing the overlap structure of clone maps. Priorknowledge of consecutive positives shrinks the sizeof the sample space of positives, which hence reduces
the size of pooling design. This observation is consistent with the original motivation of Balding andTorney [4] and Colbourn [9].
Additionally, ILB(|Ω(dd)|, f, m, k) can be used to
predict the influence of f and k on the positive
detectability of random k-set designs. To show this,we re-examine one of the screening problems studied
by Knill et al. [27] and Uehara and Jimbo [37]. In
the problem, the DNA library consists of 1298 clones
without linear orders, among
there are four� which�
positives, and thus |Ω(dd)| = 1298.
Assigningtwo set4of values to f , we fix f1 and f2 as shown in TABLE 1
and TABLE 2, respectively.

d=4

0.5

n=12980.1

For convenience, we denote by ILB(|Ω(dd)|, f, m, k)
the information-theoretic lower bound of (XII). By
letting EA [perr (A)] = 0, (XII) derives a necessary
condition on m and k.

bound of |Ω(dd)|. This can be seen as follows. To begin
with, it is obvious that there are d1 − 2d2 maximum
consecutive subset with a single positive and d2 maximum consecutive subset with 2 consecutive positives,
and therefore we can treat each single positive asa red ball and each pair of consecutive positives asa blue ball. To calculate |Ω(dd)|, it is equivalent tocounting how many ways these d1 − d2 balls with
distinguishable colors can be put into n − (d1 − d2 ) + 1
urns such that any two ball have at least an urn
between. This constraint implies that each ball will
occupy at most 3 urns. Thus, to place the ith ball, atleast (n−(d1 −d2 )+1−3(i−1)) urns are left to �choose�
2from, for i = 1, 2, . . . , d1 − d2 . Noticing that d1d−d2
is the number of ways to color these balls without
repetition, the desired lower bound is obtained.
With |Ω(dd)| and f , we give the algorithm called RPD(random pooling designer) for exhaustive search of
the least value of m and some value of k that satisfy
the Necessity Theorem. md,f and kd,f denote by theleast value of m and by the value of k with respect

10

to |Ω(dd)| and f , respectively. With the values of md ,fand kd ,f , various strategies can be used to choose theproper values of m and k for θ̂ and f .
Random Pooling Designer
Input: |Ω(dd)| and f
Output: md ,f and kd ,f
md ,f ← 1
kd ,f ← 1
temp ← ILB(|Ω(dd)|, f, md ,f , kd ,f )while temp > 0 do
md ,f ← md ,f + 1
kd ,f ← 1
while temp > ILB(|Ω(dd)|, f, md ,f , kd ,f ) and k ≤m dotemp ← ILB(|Ω(dd)|, f, md ,f , kd ,f )kd ,f ← kd ,f + 1
end while
end while5

S IMULATIONThis section shows the performance of MMCPD withrandom k-set designs chosen by RPD.
5.1 Simulation Method
The simulations are performed as follows.
1) Letting n = 1298 be the size of C, the priorknowledge of consecutive positives is set tod1 = 6, d2 = 2.4, d3 = 0.01 and d4 = 0.001.
2) The likelihoods of experimental error f are set
to f1 given in TABLE 1.
3) Find an approximate prior probability distribution Pθ that incorporates the prior knowledge by
solving (V).
4) Estimate the distribution of positive patterns ofPθ .
5) Choose proper positive patterns, and compute|Ω(dd)|, md ,f and kd ,f , for each chosen d .
6) Choose mθ,f and kθ,f such that mθ,f is the least
value with ILB(|Ω(dd)|, f1 , mθ,f , kθ,f ) < 0, for all
the chosen positive patterns.
7) The positive set P is given in two ways.
• Fix {240, 241, 890, 891, 1001, 1002} as the positive set.
• Randomly choose a positive set approximately according to Pθ .
8) Generate a random k-set design A with someproper values of k and m.
9) Compute the number of positives in each pool.
10) Based on 11), determine the pooling results r
randomly according to f1 .
11) Implement MMCPD to decode the corrupted
pooling results r. The posterior probability of
being a positive is estimated for each clone.
12) Clones are sorted in a decreasing order according to their posterior probabilities.

5.2

Preprocess of Pooling Procdedure

Solving (VII) with Descombes et al. [11]’s method,
we find that θ̂ = (6.96, −7.65, 4.955, 2.01) is a suitableapproximation for the the desired prior distribution.
However, the estimation based on MCMC simulationis costly and hence a good initial θ will accelerate
the convergence of Descombes et al [11]’s method.
In implementation, to find a proper initial θ involves
some guesswork and it can be done in a bisection way
by using trial and failure method.
Next, we take 10000 samples for estimating the distribution of positive patterns of Pθ̂ and list the positive
patterns with sample mean above 0.0025 in TABLE 3.
It also lists the corresponding md,f s and kd,f s withrespect to the positive patterns of our considerationand f1 . Particularly, d = (0, 0, 0, 0) is trivial and weremove it from our consideration. The total positive
patterns of our consideration roughly account for 95%
of the positive patterns of the samples.
TABLE 3: Estimated Distribution of d with respect toPθ̂ and Values of md ,f and kd ,f with respect to f1
dSample Mean

md,f

kd,f

(0, 0, 0, 0)(1, 0, 0, 0)(2, 0, 0, 0)(2, 1, 0, 0)(3, 0, 0, 0)(3, 1, 0, 0)(4, 1, 0, 0)(4, 2, 0, 0)(5, 1, 0, 0)(5, 2, 0, 0)(6, 1, 0, 0)(6, 2, 0, 0)(6, 3, 0, 0)(7, 2, 0, 0)(7, 3, 0, 0)(8, 2, 0, 0)(8, 3, 0, 0)(8, 4, 0, 0)(9, 3, 0, 0)(9, 4, 0, 0)(10, 3, 0, 0)(10, 4, 0, 0)(10, 5, 0, 0)(11, 4, 0, 0)(11, 5, 0, 0)(12, 4, 0, 0)(12, 5, 0, 0)(12, 6, 0, 0)(13, 5, 0, 0)(13, 6, 0, 0)(14, 6, 0, 0)(14, 7, 0, 0)(15, 7, 0, 0)
Others

0.02860.03070.01920.065
0.00750.07980.04700.07910.0191
0.09250.006
0.0617
0.06180.02320.07170.00710.04590.03850.01890.04460.00480.02720.01810.01180.02170.00410.01280.00590.00400.00870.00480.00280.00380.0216

17
33
17
48
33
48
31
63
46
96
61
44
75
58
88
72
55
85
69
98
82
65
95
78
108
91
74
104
86
99
82
95
-

7
9
4
10
6
8
5
8
6
9
7
5
7
6
8
6
4
7
5
7
6
4
6
5
6
5
4
6
5
5
4
4
-

11

TABLE 4: Positive Detectability of MMCPD for Fixed
Positive Set|P ∩ CP |
6

Times
98

5

39

4

149

3

37

2

112

1

23

0

42

During the simulations, we observed that MMCPD
detects the underlying true positives in a pairwise
way, which can also be seen from the TABLE 4.
Successfully decoding one underlying true positive
will remarkably improve the performance of detecting
the consecutive one. Besides, we can also see thatthe lower bound ILB(|Ω(dd)|, f1 , md ,f , kd ,f ), though
underestimates the underlying true minimal numberof pools required, is nontrivial and useful.
5.4

Simulation 2: Random Positive Set

By using the TABLE 3, we can verify that mθ̂,f = 108
is the least value with ILB(|Ω(dd)|, f1 , mθ̂,f , kθ̂,f ) < 0,for all the chosen positive patterns, where kθ̂,f = 6.
1000 simulations are implemented to show MMCPD’s positive detectability. In each simulation, anonempty positive set is first randomly generated
approximately according to Pθ̂ and then a random 6set design of size 108 is independently constructed.
The decoding procedure of random positive set is

0.8
0.6
0.4

Success Ratio Suc(P)0.2

Of all the chosen positive patterns, d = (6, 3, 0, 0)
attracts our attention because it has relatively high frequency, contains relatively large number of positives,
but requires few pools with small value of k. By fixing
the positive set P = {240, 241, 890, 891, 1001, 1002},
500 simulations are implemented to show MMCPD’s
positive detectability with random 5-set designs ofsize 44.
In each simulation, beginning with the state drawn
6from i.i.d Bernoulli trials with parameter q = 1298,
the warmup period includes 5000 steps. After the
warmup period, another 15000 steps are run, and thestates obtained in every 3 steps are used as samples
to estimate the posterior probability of being positivefor each clone, that is, approximately the proportionof the obtained states including the given clone. We
subsample the Markov chain in hope of weakening
potential autocorrelations among the samples. Denoteby CP the set of candidate positives. MMCPD outputs
CP which consists of the six clones with the highest mean posterior probabilities. TABLE 4 shows thenumber of times among 500 simulations that |P ∩ CP |positives can be identified by using MMCPD.

1.0

Simulation 1: Fixed Positive Set

0.0

5.3

0

200

400

600

800

1000

Ordered Simulation Index

Fig. 2: Detectability of MMCPD for Random Positiveset

the same with that of the fixed positive set except
for the choice of the set of candidate positives. To
evaluate the detectability of MMCPD for randomly
generated positive set P , we introduce the success|ratio Suc(P ) � |P ∩CP|P | , where CP is the candidate
positive set of |P | clones which have the highest
mean posterior probabilities. Obviously, the complete
identifiability occurs when Suc(P ) = 1. This measurement is nevertheless strict. Fig. 2 shows the simulation
results.
In the simulations, MMCPD showed a stable andreliable performance with random designs. Its novel
detectability is due to the prior knowledge of consecutive positives and also due to the proper choice ofvalues of m and k suggested by RPD. There are 593
instances where MMCPD completely identified all the
randomly generated positives. However, if too many
positives or/and errors occur, MMCPD tends to fail.
Denote by EF N and EF P the number of false negative
errors and that of false positive errors, respectively.
TABLE 5 lists all the instances with the success ratio
less than 0.4.
TABLE 5: Instances Where MMCPD Performed
Badly
Index

d

EF N

EF P

Suc(P )

1

(1, 0, 0, 0)

0

21

0

2

(2, 0, 0, 0)

4

17

0

3

(1, 0, 0, 0)

0

14

0

4

(15, 6, 0, 0)4

70.133333

5(21, 10, 0, 0)7

60.285714

6(3, 0, 0, 0)1

130.333333

7(3, 0, 0, 0)2

180.333333

8(12, 5, 0, 0)1

130.333333

9(11, 4, 0, 0)6

70.363636

10(16, 7, 0, 0)

4

4

0.375

11

(8, 3, 0, 0)

3

17

0.375

12

Typically, these instances are hard. If too many
positives or/and experimental errors occur, or if thepositive pattern of consecutive positives deviates far
from our prior knowledge, or a concurrence of someof these difficulties, will weaken the positive detectability of MMCPD. To overcome this, we mayeither enlarge the set of candidate positives or use
a pooling design with more pools, or both.

6

In this paper, we study the problem of screening
clone maps in the presence of experimental errors. In
conclusion, as was foreseen in Balding and Torney [4]the overlap structure of clone maps facilitates efficientpooling experiments. Our method also allows non
unique-screening as long as the prior knowledge canappropriately be interpreted.

A PPENDIX A
P ROOF OF L EMMA 1��
EA HA (rh ) can be written into
� �
�
3
�
Pr(A) −PrA (rh = j) log PrA (rh = j) .

Notice that −x log x defined on R+ is a concave
function. Hence, by using Jensen’s
inequality
we can��obtain an upper bound of EA HA (rh ) , that is,
EA HA (rh )=�

j=0 A∈Am,k

≤

�

j=0

A∈Am,k· log� �Pr(A)PrA (rh = j)A∈Am,k�AΩ(dd) is fixed and known, due to the fact that X −→f
AX −→ r forms a Markov chain we have�
PrA (rh = j) =PrA (X)PrA (rh = j|X).
d)
X∈Ω(d

Due to Assumptions 3 and 4, we have
PrA (rh = j|X)

=
=

1
�
i=0
1
�
i=0
1
�
i=0

PrA (rh = j, (AX)h = i|X)
Pr(rh = j|(AX)h = i)Pr((AX)h = i|X)f (j, i)Pr((AX)h = i|X).

Pr(A)A∈Am,k

�

�

d)
X∈Ω(d

Pr(X)

�

PrA (X)PrA (rh = j|X) (A. 2)
Pr(A)PrA (rh = j|X).

A∈Am,kd)X∈Ω(dSince for any fixed X with |dd| positives, Pr((AX)h =i|X) = 1 if pooling design A contains at least one
positive in the hth pool, otherwise 0. When the hthpool is considered, if X is fixed and A is random kset design constructed by independent column generation of uniform random k-sets and hence A ischosen independently and uniformly from Am,k , thenk |dPr((AX)h = 0|X) = 1 with probability (1 − m) d| .
Therefore, we have�Pr(A)PrA (rh = j|X)

A∈Am,k

=

1
�
i=0

1
�

f (j, i)�

Pr(A)Pr((AX)h = i|X)

A∈Am,k�k=f (j, i) 1 −mi=0

�|dd|(1−i) �

�

k
1− 1−
m

�|dd| �i.

(A. 3)A PPENDIX B
P ROOF OF L EMMA 2

�Hence, to end our proof it suffices to calculateA∈Am,k Pr(A)PrA (rh = j). Since the sample space

=

�

=

Consequently, we finish the proof by combining theresults of (A. 1), (A. 2) and (A. 3).

− Pr(A)PrA (rh = j) log PrA (rh = j)

� �
�
3
�
−
Pr(A)PrA (rh = j)

Pr(A)PrA (rh = j)A∈Am,kj=0A∈Am,k3
��

=

C ONCLUSION�

Since PrA (X) = Pr(X) is assumed in Assumption 7,
we further have(A. 1)

We assume |dd|k ≤ m. For a fixed XP with |dd| positivesd||dand a fixed t, let AXP ,m,k (t) be the subset of all
the k-set designs of size m satisfying |AXP | = t.
Obviously, the value of t can range from k to |dd|k,and it is determined by the columns corresponding to
the positive set P . Since the columns of random k-set
designs are randomly and independently generated,d||d
it follows that PXP ,m,k (t) only depends on |dd|, m andk. Therefore, any XP with |dd| positives can be used
d|
d|
|d
|d
to derive Pm,k (t). Let uXP ,m,k (t) be the number ofways of choosing columns corresponding to XP suchthat any k-set design with those columns belongs tod||d
AXP ,m,k (t). Due to construction method of random
k-set designs, we have
d|
d|
|d
|d
|AXP ,m,k (t)|uXP ,m,k (t)d||d� m �nPm,k (t) ==�m�d .
kk

Hence, to end our proof it suffices to give and||diterative expression of uXP ,m,k (t). Without loss ofgenerality and for the sake of simplicity, we express

13
d|
|d

d−1|
|d

uXP ,m,k (t) by uXP � ,m,k (t� )s, for some P � ⊂ P with|dd| − 1 positives.
Let v1 , . . . , v|dd|−1 be the columns of a k-set design A corresponding the positive set P � and v|dd|
be the column corresponding to P \ P � . Notice that,for w = 0, 1, . . . , k, if there are t − k + w positive
�|dd|−1
coordinates in the vector
j=1 vj , then we have�t−k+w��m−(t−k+w)�
ways to choose v|dd| such thatk−w
�|ddw|−1v+v
hast
positive
coordinates. This implies
j
d|
|d
j=1
that
�
�
k
�
t−k+w
d|
|d
uXP ,m,k (t) =δ(t − k + w, k)
w
w=0
�
�
· m − (tk −− wk + w) u|dXd|−1(t − k + w).
P � ,m,k� �|dd|Dividing by mon both sides, the proof ends.
k

A PPENDIX CP ROOF OF L EMMA 3
By the definition of conditional entropy and Assumption 7, we have��
EA HA (r|X)
�
�
=
Pr(A)
PrA (X)HA (r|X)
A∈Am,k

=

�

d)
X∈Ω(d

�

Pr(X)

Pr(A)

A∈Am,k

d)
X∈Ω(d

�

�

�

r∈{0,1,2,3}m−PrA (r|X)· log PrA (r|X) .

(C. 1)

Given any X ∈ Ω(dd), we have the partition
Am,k =

d|k
|d

�

t=k

d|
|d

AX,m,k (t).

Moreover, recall that PrA (r|X) =m�

f (rh , (AX)h )h=1

due to Assumptions 3 and 4. Therefore, this impliesd||dwhen X ∈ Ω(dd) is fixed, for any given A ∈ AX,m,k (t),HA (r|X) =m�

H(rh |(AX)h )

h=1
3
�

=t

j=0

−f (j, 1) log f (j, 1)+ (m − t)3
�
j=0(C. 2)−f (j, 0) log f (j, 0)= tHf (1) + (m − t)Hf (0).
Using the result of (C. 2) and Lemma 3, it follows that

(C.1) =

�

Pr(X)

d|k
|d
�

�

Pr(A)

d|
t=k A∈A|d

d)
X∈Ω(d

�

X,m,k (t)· tHf (1) + (m − t)Hf (0)
=

�

Pr(X)

d|k
|d
�

d|
|d
Pm,k (t)

t=k

d)
X∈Ω(d

+ (m − t)Hf (0)
=

d|k
|d
�

d|
|d
Pm,k (t)

t=k

�

�

�

�

tHf (1)

�tHf (1) + (m − t)Hf (0) .

ACKNOWLEDGMENTS
The authors deeply thank Professor Maiko Shigeno
for enlightening discussions and helpful comments.

R EFERENCES[1] G. Atia and V. Saligrama, ”Boolean compressed sensing andnoisy group testing,” CoRR, vol. abs/0907.1061, 2009.
[2] E. Barillot, B. Lacroix, and D. Cohen, ”Theoretical Analysis of
Library Screening Using an N-Dimensional Pooling Strategy,”
Nucleic Acids Research, vol. 19, no. 22, pp. 6241-6247, 1991.
[3] T. Berger, J.W. Mandell, and P. Subrahmanya, ”Maximally Efficient Two-Stage Screening,” Biometrics, vol. 56, no. 3, pp. 833840, 2000.
[4] D.J. Balding and D.C.Torney, ”The design of pooling experiments for screening a clone map,” Fungal Genet. Biol., 21,
pp.302-307, 1997.
[5] W.J. Bruno, E. Knill, D.J. Balding, D.C. Bruce, N.A. Doggett,
W.W. Sawhill, R.L. Stallings, C.C. Whittaker, and D.C. Torney,”Efficient Pooling Designs for Library Screening,” Genomics,
vol. 26, no. 1, pp. 21-30, 1995.
[6] W.J. Bruno, F. Sun, D.C. Torney, ”Optimizing nonadaptive
group tests for objects with heterogenous priors,” SIAM J. Appl.
Math., 58, pp. 1043-1059, 1998.
[7] C. L. Chan, P. H. Che, S. Jaggi, and V. Saligrama, ”Non-adaptive
probabilistic group testing with noisy measurements: Nearoptimal bounds with efficient algorithms,” in Communication,
Control, and Computing (Allerton), 2011 49th Annual Allerton
Conference on, Oct. 2010.
[8] C. L. Chan, S. Jaggi, V. Saligrama, and S. Agnihotri, ”Nonadaptive group testing: Explicit bounds and novel algorithms,”
CoRR, vol. abs/1202.0206, 2012.
[9] C.J. Colbourn, ”Group testing for consecutive positives,” Ann.
Combinatorics, 3, pp. 37-41, 1999.
[10] T. Cover and J. Thomas, Elements of Information Theory. John
Wiley and Sons, 1991.
[11] X. Descombes, R. Morris, J. Zerubia, and M. Berthod, ”Maximum likelihood estimation of Markov random field parameters
using markov chain Monte Carlo algorithms,” Proc. Int. Workshop EMMCVPR’97: Energy Minimization Methods Computer
Vision Pattern Recognition, vol.1223, pp.133-148, 1997.
[12] D.Z. Du and F.K. Hwang, Combinatorial Group Testing and Its
Application. World Scientific, 2000.
[13] D.Z. Du and F.K. Hwang, Pooling Designs and Nonadaptive
Group Testing: Important Tools for DNA Sequencing. World Scientific, 2006.
[14] Y.X. Cheng and D.Z. Du, ”New Constructions of One- andTwo-Stage Pooling Designs,” J. Computational Biology, vol. 15,
no. 2, pp. 195-205, 2008.
[15] R. Dorfman, ”The Detection of Defective Members of Large
Populations,” Ann. Math. Statistics, vol.14, No.4, pp. 436-440,
1943.

14

[16] A.G. D’yachkov, F.K. Hwang, A.J. Macula, P.A. Vilenkin, andC.Weng, ”A Construction of Pooling Designs with Some Happy
Surprises,” J. Computational Biology, vol. 12, no. 8, pp. 11291136, 2005.
[17] R.G. Gallager, ”Low-Density Parity-Check Codes,” IRE Trans.
Information Theory, vol. 8, pp. 21-28, 1962.
[18] G. Ge, Y. Miao, and X. Zhang, ”On Block Sequences of
Steiner Quadruple Systems with Error Correcting Consecutive
Unions,” SIAM J. Discrete Math., pp. 940-958, 2009.
[19] S. Geman and D. Geman, ”Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images,” IEEE Trans.
Pattern Anal. Machine Intell., vol. PAMl-6, pp. 721-741, 1984.
[20] C.J. Geyer and E. A. Thompson, ”Constrained Monte Carlomaximum likelihood for dependent data (with discussion),” J.
Roy. Statist. Soc., vol.54, pp. 657-699, 1992.
[21] E.T. Jaynes, ”Information theory and statistical mechanics,”
Phys. Rev. 106, pp. 620-630, 1957.
[22] E.T. Jaynes, ”Information theory and statistical mechanics, II,”
Phys. Rev. 108, pp. 171-190, 1957.
[23] E.T. Jaynes, ”Prior Probabilities,” IEEE Trans on Systems Science and Cybernetics, vol 4, no. 3, pp. 227-241, 1963.
[24] J.S. Juan, G.J. Chang, ”Adaptive group testing for consecutive
positives,” Discrete Math., vol. 308, 7, pp. 1124-1129, 2008.
[25] R. Kindermann and J.L Snell, Markov Random Fields and Their
Applications. American Mathematical Society, 1980.
[26] E. Knill, ”Lower bounds for identifying subset members withsubset queries,” Proc. 6th ACM-SIAM sump. Discrete Algorithms, pp 369-377, 1995.
[27] E. Knill, A. Schliep, and D.C. Torney, ”Interpretation of PoolingExperiments Using the Markov Chain Monte Carlo Method,”
J. Computational Biology, vol. 3, no. 3, pp. 395-406, 1996.
[28] D.J.C. MacKay and R.M. Neal, ”Near Shannon Limit Performance of Low Density Parity Check Codes,” IEE Electronics
Letters, vol. 32, no. 18, pp. 1645-1655, 1996.
[29] A.J. Macula, ”Probabilistic Nonadaptive Group Testing in the
Presence of Errors and DNA Library Screening,” Ann. Combinatorics, vol. 3, no. 1, pp. 61-69, 1999.
[30] M. Mézard and C. Toninelli, ”Group Testing with Random
Pools: Optimal Two-Stage Algorithms,” IEEE Trans. Info. Th.
57, 1736-1745, 2011.
[31] M. Müller and M. Jimbo, ”Consecutive positive detectable
matrices and group testing for consecutive positives,” Discrete
Math., vol. 279, pp. 369-381, 2004.
[32] M. Müller and M. Jimbo, ”Cyclic sequences of k-subset with
distinct consecutive unions,” Discrete Math., vol. 308, pp. 457464, 2008.
[33] H.Q. Ngo and D.Z. Du, ”A Survey on Combinatorial Group
Testing Algorithms with Applications to DNA Library Screening,” Discrete Math. Problems with Medical Applications, DIMACS Ser. Discrete Math. and Theoretical Computer Science,
vol. 55, pp. 171-182, 2000.
[34] J. Pearl, Probabilistic Reasoning in Intelligent Systems: Networks
of Plausible Inference. Morgan Kaufmann, 1997.
[35] G.O. Roberts and A.F.M. Smith, ”Simple conditions for theconvergence of the Gibbs sampler and Metropolis-Hastings
algorithms,” Stochastic Process. Appl. vol. 49, pp. 207-216, 1994.
[36] L. Tierney, ”Markov Chains for exploring posterior distributions, ” Ann. Statistics, vol. 22, No. 4, pp. 1701-1728, 1994.
[37] H. Uehara and M. Jimbo, ”A positive detecting code and its
decoding algorithm for DNA library screening,” IEEE/ACM
Trans Computational Biology and Bioinformatics, vol. 6, no. 4,
pp. 652-666, 2009.
[38] P. Sham, J.S. Bader, I. Craig, M. O’Donovan, and M. Owen,”DNA Pooling: A Tool for Large-Scale Association Studies,”
Nature Rev. Genetics, vol. 3, no. 11, pp. 862-871, Nov. 2002.
[39] C.E. Shannon, ”A Mathematical Theory of Communication,”
Bell System Technical Journal, 27, pp. 379-423 & 623−656, July
& October, 1948.
[40] W. Wang, M.J. Wainwright, K. Ramchandran, ”InformationTheoretic Limits on Sparse Signal Recovery: Dense Versus
Sparse Measurement Matrices,” arXiv:0806.0604, 2008.
[41] M.J. Wainwright, ”Information-theoretic bounds on sparsity
recovery in the high-dimensional and noisy setting,” IEEE
Trans. Inf. Theory, vol. 55, pp.5728-5741, 2009.