Oguz Kaya

Bora Uçar

High-performance parallel algorithms for the Tucker decomposition of higher order sparse tensors

Keywords:

We investigate an efficient parallelization of a class of algorithms for the well-known Tucker decomposition of general N -dimensional sparse tensors. The targeted algorithms are iterative and use the alternating least squares method. At each iteration, for each dimension of an N -dimensional input tensor, the following operations are performed: (i) the tensor is multiplied with (N -1) matrices (TTM step); (ii) the product is then converted to a matrix; and (iii) a few leading left singular vectors of the resulting matrix are computed (SVD step) to update one of the matrices for the next TTM step. We propose an efficient parallelization of these algorithms for current supercomputers comprised of compute nodes, where each node is a multi-core system. We reformulate the computation of N successive TTM-steps to increase the reuse of intermediate computation, which is of interest on its own. We discuss a set of preprocessing steps which takes all computational decisions out of the main iteration of the algorithm and provide an intuitive row-wise shared-memory parallelism for the TTM and SVD steps. We consider a coarse and a fine grain computational scheme, investigate their data dependencies, and identify efficient communication schemes. We demonstrate how the computation of singular vectors in the SVD step can be carried out efficiently following the TTM step. Finally, we develop a hybrid MPI-OpenMP based implementation of the overall algorithm and report speedup results on up to 2048 cores.

Algorithmes parallèles haute performance pour la décomposition Tucker des tenseurs creux

Résumé : Nous nous intéressons à la parallélisation efficace d'une classe d'algorithmes pour la décomposition Tucker de tenseurs creux en dimension N . Les algorithmes ciblés sont itératifs et utilisent la méthode alternée des moindres carrés. A chaque itération, pour chaque dimension d'un tenseur de dimension N, les opérations suivantes sont effectuées: (i) le tenseur est multiplié par N-1 matrices (étape TTM); (ii) le produit est converti en une matrice; et (iii) des vecteurs gauches singuliers de la matrice résultante sont calculés (étape SVD) afin de mettre à jour l'une des matrices pour la prochaine étape TTM. Nous proposons une parallélisation efficace de ces algorithmes pour les supercalculateurs dans lesquels chaque noeud est un système multi-coeurs. Nous reformulons le calcul de N étapes TTM successives afin d'améliorer la réutilisation de calculs intermédiaires. Nous discutons un ensemble d'étapes de pré-traitement qui sortent les décisions de calcul hors de l'itération principale de l'algorithme, et nous fournissons une parallélisation intuitive pour les étapes TTM et SVD. Nous considérons un schéma de calcul à grain fin et à gros grain, nous étudions les dépendances de données, et nous identifions des schémas de communication efficaces. Nous démontrons comment le calcul de vecteurs singuliers dans l'étape SVD peut être menée efficacement suite à l'étape TTM. Finalement, nous développons une implémentation hybride MPI-OpenMP de l'algorithme global et nous obtenons les facteurs d'accélération avec jusqu'à 2048 coeurs.

Mots-clés : tenseurs creux, factorizations des tenseurs creux d'ordre large

Introduction

Tensors or multi-dimensional arrays are used to represent data with higher dimensionality in many applications. Some of the most popular applications include analysis of Web graphs [START_REF] Kolda | The TOPHITS model for higher-order web link analysis[END_REF], knowledge bases [START_REF] Carlson | Toward an architecture for never-ending language learning[END_REF], product reviews at online stores [START_REF] Bennett | The netflix prize[END_REF], chemometrics [START_REF] Appellof | Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents[END_REF], signal processing [START_REF] Lathauwer | From matrix to tensor: Multilinear algebra and signal processing[END_REF], computer vision [START_REF] Vasilescu | Multilinear analysis of image ensembles: Tensorfaces[END_REF], forensics [START_REF] Maruhashi | MultiAspectForensics: Pattern mining on largescale heterogeneous networks with tensor analysis[END_REF] and more. Tensor decomposition algorithms are used to find latent relations in the tensors. There are two prominent tensor decomposition formulations. CANDE-COMP/PARAFAC (CP) formulation approximates a tensor as a sum of the product of rank-one tensors. Tucker formulation approximates a tensor with a core tensor multiplied by a matrix along each dimension or mode-see Fig. 1 for a simplistic view. Both of these formulations have uses in applications; in particular the CP formulation is deemed useful respectively for understanding latent factors, and the Tucker formulation is deemed more appropriate for compressions and identifying relation among the factors [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF]. In this work, we investigate efficient computation of Tucker decomposition of sparse high dimensional tensors in common distributed memory systems.

There are variants of CP and Tucker decompositions, and different algorithms to compute them [START_REF] Comon | Tensors: A brief introduction[END_REF][START_REF] Kolda | Tensor decompositions and applications[END_REF]. The most common algorithms for both decompositions and their variants are based on the alternating least squares (ALS) method. The algorithms of this type are iterative, where the computational core of an iteration is a special operation among an N -mode tensor and N matrices. The key operation in the ALS-based CP decomposition (CP-ALS) case is called the matricized tensor times Khatri-Rao product (MTTKRP); we refer the reader to other resources for details [START_REF] Kolda | Tensor decompositions and applications[END_REF]. The key operation in the ALS-based Tucker decomposition (Tucker-ALS) method is called the tensor times matrix-chain (TTMc) product. These two operations pose similar challenges for efficient computations. But there are distinct challenges and opportunities for efficiency as well.

The tensor times matrix-chain (TTMc) product in the Tucker-ALS method is performed for all modes of the N -mode input tensor at every ALS iteration. TTMc for a mode n involves tensor times matrix (TTM) products with N -1 different matrices, each of which is associated with one of the modes other than n. The TTM product can be considered as a higher dimensional variant of the matrix-vector product operation (Section 2 explains the TTM operation in detail). Techniques for efficiency of a single TTM are therefore akin to those in the matrix-vector case but require increased effort to overcome the difficulties associated with the higher dimensionality. The chain product TTMc adds another level of complexity; one needs to be careful about the memory requirements. Furthermore, one needs to take the content of an ALS iteration into account as N TTMc's will follow another.

The Tucker-ALS algorithm also computes a few singular vectors of a large, usually tall-andskinny dense matrix for each mode at every ALS iteration. This matrix arises by a logical reorganization of the result of the TTMc associated with the corresponding mode. The cost of computing the singular vectors is not negligible and hence needs to be addressed in an efficient parallel algorithm for the Tucker decomposition. We refer to the computation of the desired singular vectors as the SVD-step.

Our contributions are as follows. First, we design efficient parallel algorithms for the TTMc operation. Here our analysis points that with proper task definitions, the computational and communication requirements are as in the MTTKRP operation. Having identified this, we make use of the hypergraph models from our earlier work on CP-ALS [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF] for reducing communication and achieving load balance during each TTMc operation. Second, in order to increase the data reuse between the successive TTMc's, we propose a reformulation by generalizing exiting memory-

X ⇡ A G B C Figure 1 -Tucker decomposition of a 3rd mode tensor X ∈ R I 1 ×I 2 ×I 3 as a core tensor G ∈ R R 1 ×R 2 ×R 3 multiplied by matrices A ∈ R R 1 ×I 1 , B ∈ R R 2 ×I 2 and C ∈ R R 3 ×I 3 in different modes.
In CP-decomposition, G is a diagonal tensor having the same size along each dimension, A, B and C have the same number of rows.

efficient implementations of TTMc's [6,[START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF]. This is of independent interest and can be used with dense or sparse tensors. Third, we propose efficient realization of the computation of the singular vectors by making use of existing libraries PETSc [START_REF] Balay | PETSc Web page[END_REF] and SLEPc [START_REF] Roman | SLEPc users manual[END_REF]. We carefully designed this step so that the communication requirements in parallel iterative algorithms used for computing the singular vectors are reduced and the load balance is achieved by making use of the data decomposition of the TTMc's step. Fourth, we put every thing together and propose an OpenMP-MPI hybrid and complete implementation of Tucker-ALS algorithm and present speed-ups on a high-end parallel system where up to 742 speed-up on 2048 cores are reported on real world data. The organization is as follows. We give background on the basic tensor operations in the next section. Then, in Section 3, we discuss the proposed parallel Tucker-ALS algorithm. Here, we first discuss the reorganization of TTMc's for efficient data reuse. We then discuss parallel TTMc's and SVD-steps. We then give a brief summary of related recent work in Section 4. Experimental results on a parallel system with compute nodes containing 32 cores are presented in Section 5. We conclude the paper with a few remarks in Section 6.

Background

We use bold, upper case Roman letters for matrices, as in A. The matrix elements are shown with the corresponding lowercase letters, as in a i,j . Matlab notation is used to refer to the entire rows and columns of a matrix, e.g., A(i, :) and A(:, j) refer to the ith row and jth column of A respectively.

We use calligraphic font to refer to tensors, e.g., X. The order of a tensor is the number of its dimensions or modes, which we denote with N . For the sake of simplicity of the notation and the discussion, we sometimes discuss the case N = 3, even though our algorithms and implementations have no such restriction. We explicitly generalize the discussion to general order-N tensors whenever we find necessary. As in matrices, an element of a tensor is denoted by a lowercase letter and subscripts corresponding to the indices of the element, e.g., the element (i, j, k) of a third-order tensor is x i,j,k . A fiber in a tensor is defined by fixing every index but one, e.g., if X is a third-order tensor, X :,j,k is a mode-1 fiber and X i,j,: is a mode-3 fiber.

We reproduce the following definitions from Kolda and Bader's survey [START_REF] Kolda | Tensor decompositions and applications[END_REF]. The n-mode matricization of a tensor X is denoted by X (n) and refers to the reordering of X's elements into a matrix by arranging the nth mode fibers as the columns of X (n) . For example, take

X ∈ R I 1 ו••×I N .
Then X (1) denotes the mode-1 matricization of X where the rows of X (1) correspond to the first mode of X and the columns correspond to the remaining modes. The tensor element x i 1 ,...,i N corresponds to the element i 1 , 1 + N j=2 (i j -1) j-1 k=1 I k of X (1) . Specifically, each column of the matrix X (1) is a mode-1 fiber of the tensor X.

The n-mode product of a tensor X ∈ R I 1 ו••×I N with a matrix U ∈ R J×In is denoted by X × n U. This is also referred to as tensor times matrix (TTM) product. The result Y is a tensor of size

I 1 × • • • × I n-1 × J × I n+1 × • • • × I N . A particular element of Y is given by y i 1 ,...,i n-1 ,j,i n+1 ,...,i N = In in=1 x i 1 ,i 2 ,...,i N u jin .
A tensor can be multiplied by a set of matrices along a given set S of modes. We use the notation TTMc(X, S, {U n : for n ∈ S}) to refer to the tensor n-mode product of X with matrices U n for n ∈ S. We use TTMc(S) for clarity, as the tensor X and the matrices U n 's will be clear from the context.

The Tucker decomposition expresses a given tensor X ∈ R I 1 ו••×I N as a core tensor G multiplied by a factor matrix U n of size I n × R n at each mode n. Here, R 1 , . . . , R N are the requested rank of the decomposition for each mode. Formally, the Tucker decomposition

[[G; U 1 , . . . , U N ]] writes X as G × 1 U 1 × 2 • • • × N U N . For example, if X ∈ R I 1 ×I 2 ×I 3 , then its Tucker decomposition [[G; A, B, C]] writes x i,j,k ≈ I 1 p=1 I 2 q=1 I 3
r=1 g pqr a ip b jq c kr . A well-known algorithm for computing the Tucker decomposition is based on the alternating least squares method which is given in Algorithm 1, which is also called higher-order orthogonal iteration (HOOI) [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. In this algorithm, the factor matrices are initialized first. This initialization could be done randomly or using the higher-order SVD [START_REF] Lathauwer | A multilinear singular value decomposition[END_REF]. Then, the alternating least squares method is applied in the "repeat-until" loop. Here, for each mode n, the tensor matrix-chain multiplication TTMc({1, . . . , N } \ {n}) at Line 1 is computed. This yields a tensor of size

R 1 × R 2 × • • • R n-1 × I n × R n+1 × • • • × R N , which is then matricized along the nth mode yielding Y (n) ∈ R In×Π i =n R i ,
and the left R n leading singular vectors of Y (n) are computed at Line 2 and used as the columns of U n .

The algorithms for computing the singular vectors are well documented [START_REF] Golub | Matrix Computations[END_REF]Ch. 8.6]. The TTMc operation at Line 1 needs attention. This is because of the fact that there are many ways to compute the product [START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF], as the TTM's can be performed in any order [START_REF] Kolda | Tensor decompositions and applications[END_REF]. One extreme, which is called the standard way, is to compute one n-mode product at a time, e.g., for a 3rd order

RR n°8801

Algorithm 1: Tucker-ALS for N -mode tensors Input : X: An N -mode tensor R 1 , . . . , R N :

The rank of the decomposition for each mode Output:

Tucker decomposition [[G; U 1 , . . . , U N ]] Initialize the matrix U n ∈ R In×Rn for n = 1, . . . , N repeat for n = 1, . . . , N do 1 Y ← X × 1 U T 1 • • • × n-1 U T n-1 × n+1 U T n+1 • • • × N U T N 2 U n ← R n leading left singular vectors of Y (n) until no improvement or maximum iterations reached G ← X × 1 U T 1 • • • × N U T N return [[G; U 1 , . . . , U N ]] tensor and n = 1, Y = ((X × 2 U T 2 ) × 3 U T 3 ) . (1)
The other extreme, which is called the element-wise, is to compute the fibers of Y at a time, e.g., with the same example as above,

y :,j,k = X × 2 U 2 (j, :) × 3 U 3 (k, , :) , (2)
for j = 1, . . . I 2 and k = 1, . . . , I 3 . In the element-wise computation (2), the n-mode products are with the row-vectors of the factor matrices, where the unit size in the nth mode of the result is not indexed. Kolda and Sun also discuss schemes for higher order tensors, where some of the modes are treated element-wise and some others in the standard way and present formulas to decide this for a given tensor. We re-write the elementary-way (2) in Algorithm 2, as it better suits for our presentation. In this algorithm, U 2 (j, :) • U 3 (k, :) refers to the outer product of the vectors U 2 (j, :) T and U 3 (k, :).

Algorithm 2: Elementary-way of multiplication in TTMc for 3rd order tensors and n = 1

foreach x i,j,k ∈ X do Y(i, :, :) ← Y(i, :, :) + x i,j,k (U 2 (j, :) • U 3 (k, :))

Parallel Tucker Factorization

Here we propose a computational scheme to compute a set of tensor times matrix chain products efficiently. We then discuss our parallelization approach for TTMc's and for computing the singular vectors.

TTM-tree

In an iteration of the Tucker-ALS, N TTMc's are performed one after another, each one for a mode potentially having overlaps. For example, TTMc({1, . . . N }\{N }) and TTMc({1, . . . N }\{N -1})

have N -2 TTM's in common. This has been exploited before [6] to increase data reuse (see our discussion in Section 4). Here, we generalize this idea and increase the data reuse considerably. We propose a computational scheme that we call TTM-tree. This tree represents N TTMc's compactly and enable the re-use of the common TTM's. Each node in a TTM-tree corresponds to a tensor with a set of multiplied and a set of unmultiplied modes which partition {1, . . . , N }. The tensor associated with a tree node corresponds to the result of the TTMc of the input tensor X with the node's multiplied modes. A TTM-tree always has X as its root, with an empty set of multiplied modes. Starting from the root, the unmultiplied modes of the children of each node partition the node's unmultiplied modes. A TTM-tree has exactly N leaves where the n th leaf has only the mode n in its unmultiplied mode list; its result is TTMc({1, . . . , N } \ {n}) and hence the tensor Y in Line 1 of Algorithm 1. We note that the TTM-tree is a general computational scheme that can be applied to both sparse and dense Tucker computations.

We perform two types of operations on a TTM-tree while computing TTMc({1, . . . , N } \ {n}) for all n in an iteration of Tucker-ALS. In an update operation on a particular tree node, we update the TTMc result stored in the node by performing relevant TTM operations along the tree. In doing so, we first perform an update operation on the parent in case its result is not computed, or needs to be re-computed due to a recent update to a matrix U k in one of its multiplied modes k. Note that the parent might also recursively trigger another update on its parent; thereby, all results for the nodes from the leaf to the root might potentially be updated. After ensuring that the parent's result is valid, we perform a TTMc with the parent's tensor and the relevant "missing" multiplied modes of the child node to update the tensor at the child node. When U n is updated, we perform an invalidate operation. This operation frees all tensors at tree nodes whose value depends on the old value of U n , and mark those nodes as invalid. Note that these are exactly the nodes who have n in their multiplied modes.

We investigate two types of TTM-trees we call flat and binary. A flat TTM-tree consists of the root and the leaves, with no intermediate result. This corresponds to the "element-wise" product (2). This approach minimizes the memory consumption at the cost of eliminating the potential reuse of intermediate results. In a binary TTM-tree, starting from the root, every node has two children, except the leaves. Unmultiplied modes of a node is divided into left and right halves, and the children are obtained by performing TTM in these modes.

We note the following two theorems whose relatively straightforward proofs are in the accompanying technical report [START_REF] Kaya | High performance parallel algorithms for Tucker decompositions of higher order sparse tensors[END_REF].

Theorem 1. Let X be a tensor with N = 2 k modes, where k ∈ Z ≥0 . The total number of n-mode matrix multiplications in each iteration of the Tucker-ALS algorithm using binary TTM-tree is N log 2 N .

Proof. In each iteration of the Tucker-ALS, at the k th step, the binary TTM-tree updates its k th leaf, then invalidates all its nodes whose set of multiplied modes involves k. Note that by the construction of the binary TTM-tree, each ancestor of the leaf node k involves k in its set of unmultiplied nodes; hence, it does not get invalidated by the update to k, nor to any of its other descendants. Also, descendant leaves of each interior parent node are consecutive, as the unmultiplied modes of each node consists of a set of consecutive modes, which is implied by the construction procedure described above. Therefore, as the modes are updated from 1 to N within a Tucker-ALS iteration, the parent node gets updated when its first descendant leaf is updated, retains a valid result while updates are performed on all its consecutive descendant leaves, then gets invalidated and never needed for the remaining modes; because, the leaf nodes corresponding to the remaining modes are not descendants of the parent node, and they neither incur an update operation on it nor need its result. As a result, in each iteration of Tucker-ALS, each internal node is updated exactly once. Similarly, binary TTM-tree updates each leaf node once per iteration in the main Tucker-ALS loop.

The number of n-mode matrix multiplications to be performed in each mode's update in a TTM-tree equals to the number of missing unmultiplied modes with respect to its parent. Note that in the construction of the binary TTM-tree, each node partitions its unmultiplied modes into two halves, each of which constitutes the set of unmultiplied modes of its children. Therefore, each child node misses the other half of the unmultiplied modes of its parent, hence the number of missing unmultiplied modes equals to the cardinality of the childs set of unmultiplied modes.

After these two observations, we can finalize the proof. In a Tucker-ALS iteration, each node gets updated once, and performs n-mode matrix product per its missing unmultiplied mode with respect to its parent. Note that by the construction of a binary TTM-tree, the unmultiplied modes of the nodes in its each level partition the set {1, . . . , N }; therefore, the total number of n-mode matrix multiplications performed per level is N . The binary TTM-tree has k = log 2 N levels; therefore, the total number of n-mode matrix multiplications performed within a Tucker-ALS iteration is N log 2 N . Theorem 2. Let X be a tensor with N = 2 k modes, where k ∈ Z ≥0 . Then, in any iteration of the Tucker-ALS algorithm, the number of nodes with valid results in the binary TTM-tree is at most 2 log 2 N .

Proof. At the k th step of an iteration of the Tucker-ALS, after the invalidate is performed, the results of all nodes in the binary TTM-tree whose multiplied mode involves k (or equivalently, unmultiplied modes do not involve k) are invalidated and destroyed. As stated in the Theorem 1, the unmultiplied modes of the nodes at each level of the binary TTM-tree partition the set {1, . . . , N }; therefore, only one node at each level involves k in its unmultiplied modes, hence do not get invalidated. Hence, after the invalidate is performed at any step k, the number of nodes whose results remain valid is k = log 2 N . Note also that an update operation can generate at most k = log 2 N results in the binary TTM-tree; therefore, at any instant of Tucker-ALS, the total number of nodes with valid results in the binary TTM-tree cannot exceed 2log 2 N .

In contrast to the binary TTM-tree, while using the flat TTM-tree, one needs to perform TTM's for N -1 matrices in each mode, which results in a total of N (N -1) TTM's in one iteration of the Tucker-ALS algorithm. Therefore, Theorem 1 shows that the binary TTM-tree asymptotically reduces the number of TTM's, which becomes particularly effective for higher-order tensors. The disadvantage is the required memory to store the intermediate results; Theorem 2 puts a nice upper bound on the maximum number of intermediate results at any instant of the algorithm.

TTM-tree and symbolic TTM with semi-dense tensors

We use the semi-dense tensor data structure proposed by Baskaran et al. [6] for storing tensors in the TTM-tree. A semi-dense tensor consists of sparse and dense modes, namely sm and dm, which partition {1, . . . , N }. Each entry of the semi-dense tensor consists of coordinates in sparse modes stored in idx, and a value which is a dense tensor of order |dm|. A semi-dense tensor contains symbolic information such as sparse indices of elements and the set of sparse and dense modes that would remain the same through multiple TTM calls, particularly when the tensor corresponds to a node in a TTM-tree. Therefore, it is viable to precompute and store all such symbolic information before the main Tucker-ALS iteration for efficiency. For this reason, we introduce Algorithm 3, which computes all the symbolic information needed for the numerical TTM's when executed within Tucker-ALS for many iterations.

In the symbolic TTM, from an input semi-dense tensor X p , which we call as the parent, we create the child semi-dense tensor X c holding the symbolic information for the numerical TTM of X p in the given modes. We first set the dense and sparse modes of X c in Lines 1 and 2. Note that dense modes of X c consist of the dense modes of X p and the modes, which correspond to the sparse modes of X p that become dense in X c as a result of TTM. In Line 3, the nonzero entries of X c are determined by collapsing the parent's idx array in the columns corresponding to the set of modes, which will be dense in X c , and then eliminating duplicates. Finally in Line 4, for each element i of the parent X p , we find the element j of X c to which it makes a contribution in the numerical TTM, and add it to the reduction list X c .rl(i). Thereby, we can easily perform the numerical TTM by going over the parent's elements in this update list, performing their outer product with the corresponding rows of the matrices, and sum-reducing the results in a similar fashin to the scheme provided in Algorithm 2. Note that the sparse modes of a semi-dense tensor at a node of the TTM-tree are the unmultiplied modes, whereas its dense modes are the multiplied modes. All numerical computation later is performed using the update and invalidate operations. We do not provide further details for the numerical multiplication.

Algorithm 3: Symbolic TTM for constructing semi-dense tensor X c at a child node from X p at a parent node in the TTM-tree.

Input : X p : An N -mode semi-dense tensor modes: List of sparse modes of X p to perform TTM Output: X c : The result of the TTM of X p in modes modes

1 X c .dm ← X p .dm ∪ modes 2 X c .sm ← {1, . . . , N } \ X c .dm 3 {X c .idx, X c .nz} ← collapse(X p .idx(:, X c .sm)) for i = 1, . . . , X p .nz do 4 j ← f indIndex(X c .idx, X p .idx(i, X c .sm)) 5 X c .rl(j) ← X c .rl(j) ∪ {i} return X c

Parallelization

In Algorithm 4, we provide the high-level description of the parallel Tucker-ALS algorithm at an MPI-rank p. The first step of the algorithm is to form the TTM-tree for the input tensor X p . At the beginning of each iteration, we assume that p has all the rows of U matrices that are needed to perform the TTMc operations. Having all needed data available locally, for each mode n, p performs the update operation on the leaf tensor which corresponds to rows of Y (n) owned by the MPI-rank p. In the next step, SVD is performed on the row-wise distributed matrix to obtain the new values of U n . Then, the updated rows of the U n are communicated to make the needed data available. Finally, TTM-tree invalidates and deallocates all its results that become outdated after

RR n°8801

Algorithm 4: Parallel TTM-tree based Tucker-ALS Input : X: An N -mode tensor R 1 , . . . , R N :

The rank of the decomposition for each mode Output: Tucker decomposition [[G; U 1 , . . . , U N ]] tree ← f ormT T M tree(X) repeat for n = 1, . . . , N do tree.update(tree.leaf (n))

U n ← R n leading left sing. vectors of tree.leaf (n) Send/receive the updated rows of U n tree.invalidate(n)

until no improvement or maximum iterations reached ; G ← X × 1 U T 1 • • • × N U T N return [[G; U 1 , . . . , U N ]]
U n changes. While performing updates numerically for a child X c , the ithe element is updated using the reduction lists X c .rl(i), the destinations are disjoint and hence parallel loops can be used easily to achieve node parallelism. Note that after the TTMc operation, the resulting Y tensor is matricized into

Y (n) ∈ R In×(R 1 ו••R n-1 ×R n+1 ו••×R N )
whose dimensions can be large. Therefore a parallel algorithm should be run for computing the singular vectors. SLEPc [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF] provides required algorithms. Among different alternatives [34, Ch.

4.], we chose the SVD solver based on the Lanczos bidiagonalization [START_REF] Hernández | Restarted Lanczos bidiagonalization for the SVD in SLEPc[END_REF][START_REF] Hernández | A robust and efficient parallel SVD solver based on restarted Lanczos bidiagonalization[END_REF]. This method is iterative where the computational core of each iteration is the matrix-vector (MxV) and matrix transpose-vector (MtxV) multiplications (there are some other vector operations as well) with Y (n) . By making use of PETSc's shell matrix type [START_REF] Balay | PETSc Web page[END_REF], SLEPc offers implementation choices where the user provides MxV and MtxV multiplication subroutines to be called within the SVD solvers. We parallelize the local MxV and MtxV operations at a node level by relying on the multithreaded gemv implementation in MKL. There were a few particularities to reduce the communication in the SVD solver which we discuss below at each task definition.

For the sake of the simplicity, some of the task definitions are for 3-mode tensors; however, they naturally generalize to higher dimensions.

Coarse-grain task definition

In an iteration of the Tucker-ALS, we perform two main operations for each mode n: a TTMc step to obtain a matricized tensor Y (n) , and an SVD step to obtain the matrix U n once Y (n) is ready.

In the coarse-grain task decomposition, we define computing each row i of U n as an atomic task, and hold the owner of this task responsible for computing the i th row of Y (n) . We denote this task by t n i . For n = 1, as provided in the Algorithm 2, a row i of the matricized tensor Y (n) receives a contribution x i,j,k (U 2 (j, :) • U 3 (k, :)) for each nonzero x i,j,k . Therefore, t 1 i needs all nonzeros in the tensor slice X(i, :, :), as well as the rows U 2 (j, :) and U 3 (k, :) to perform the outer product. Therefore, t 1 i may require up to |X(i, :, :)| outer products for completion (with a flat TTM-tree). Also, for each nonzero x i,j,k , t 1 i needs the data owned by t 2 j and t 3 k , which defines the data dependencies.

In our recent work [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF], we used a similar task definition in the context of CP-ALS, which we naturally represented with a hypergraph model. We can use the same model to reduce the total communication volume and to achieve load balance during the TTM computations. There is a difference though. The computational weights |X(i, :, :)|, |X(:, j, :)|, and |X(:, :, k)| for the tasks t 1 i , t 2 j , and t 3 k , respectively are exact with a flat TTM-tree, because each nonzero performs an outer product with N -1 vectors. If a binary TTM-tree is used, a nonzero may perform an outer product in all N -1 modes, or a partial outer product which gets merged with another partial outer product in an intermediate semi-dense tensor. Therefore, in this scheme, the weight becomes an upper-bound on the actual computational load.

At the end of the coarse grain TTMc, each MPI-rank has the entirety of the rows of Y (n) corresponding to its n-mode indices. The MxV multiplication y ← Y (n) x can be easily implemented by gathering all x entries to all processors as Y (n) is a dense matrix. The partition on y induces a partition on the left singular vectors of Y (n) and hence on the rows of the new U n . We aligned the partition on y with that of the rows of U n so that the owner of U n (i, :) holds y i . The MtxV x T ← y T Y (n) can be implemented by an all-to-all communication on the partial results of local y T Y (n) computations. If the MPI-ranks have almost equal number of n-mode indices, the computational load would be balanced. In our current implementation, the n-mode indices are distributed for the TTMc step among the MPI-ranks with hypergraph models to achieve load balance in terms of the number of nonzero tensor elements. This is likely to result in load balance during the MxV and MtxV operations in the SVD step.

Fine-grain task definition

We have the task u n i of computing and owning the row i of U n as in the coarse grain formulation. However, in this task definition we use a finer granularity in the TTMc operations, and define the computation of each outer product x i,j,k (U 2 (j, :) • U 3 (k, :)) as the atomic task, which we denote as z i,j,k .

For n = 1, each process which owns a task z i,j,k generates a partial result for Y [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF] . One might consider sending these partial results to the owner of the task u 1 i ; however, we avoid this communication because of the fact that each row of the matrix Y (1) has R 2 × R 3 columns, which may get quite costly even for N = 3. For N > 3, this number becomes i =1 R i , which can easily get large for high order tensors and/or ranks of approximation. Therefore, we keep the partial results in the sum-distributed form across many processors and handle this difficultly during computing singular vectors.

Modeling the tasks and their dependencies yields the same fine grain hypergraph model we proposed before [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF]. Therefore, we can reduce the total communication volume and achieve load balance during the TTM computations by making use of the results from the previous study. This hypergraph model also useful in assigning almost equal number of mode indices per MPIrank for each mode. Note that not combining the partial results of Y (n) can increase the total computational load of matrix-vector multiplications in the SVD step. This increase is equal to the λ -1 cut size metric of the corresponding hypergraph and hence is reduced. That is why the previously proposed hypergraph model is also useful for the SVD step.

At the end of the fine grain TTMc, each MPI-rank can have partial results on many nonzero elements in Y according to the distribution of the nonzeros of X. In other words, MPI-ranks share some of the rows of Y (n) where the exact value of the elements in a row can only be computed if the shared rows are added up. Since the MxVs does not really need the value of the elements of Y but Y (n) 's action on a vector, we do not accumulate those partial results on Y. Rather, we implement the MxV y ← Y (n) x for this case by reducing the partial results on y. The amount of communication on y is closely related to the communication in the fine-grain TTMc operations. Assuming a common decomposition rank of R for each mode, the total communication volume for reducing y is equal to 1/R times the total volume of communication in the TTMc, if one again aligns the partition on y with the n-mode index partition. Note that if instead one builds the exact values of Y (n) , then the same amount of communication as in the TTMc will be required. A dual communication scheme is necessary for computing the MtxVs. As in the coarse grain algorithm, the load balance in the local MxV and MtxV operations can be achieved if the MPI-ranks have equal number of rows in Y (n) . This is one of the complex partitioning problems where the total computational load can only be determined after a partition [START_REF] Kaya | On partitioning problems with complex objectives[END_REF][START_REF] Pinar | Partitioning for complex objectives[END_REF]. We do not explicitly address this problem and hope that assigning equal amount of n-mode indices will again lead to load balance. We note that the number of rows of Y (n) in which many MPI-ranks have contributions is equal to the total communication volume for reducing y, and hence again is a linear function of the TTMc's volume of communication.

Related work

We give a brief overview of mostly recent progress on efficient tensor decomposition algorithms. These can be categorized into four classes: (i) toolboxes for Matlab and similar environments [START_REF] Andersson | The N-way toolbox for MATLAB[END_REF][START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF][START_REF] Bader | Matlab tensor toolbox version 2.6[END_REF][START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF][START_REF] Papalexakis | ParCube: Sparse parallelizable CANDECOMP-PARAFAC tensor decomposition[END_REF]; (ii) implementations for shared memory systems [6,[START_REF] Baskaran | Low-overhead load-balanced scheduling for sparse tensor computations[END_REF][START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF][START_REF] Smith | SPLATT: Efficient and parallel sparse tensor-matrix multiplication[END_REF]; (iii) implementations based on MapReduce paradigm [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF][START_REF] Kang | GigaTensor: Scaling tensor analysis up by 100 times -Algorithms and discoveries[END_REF]; (iv) distributed memory systems [START_REF] Choi | DFacTo: Distributed factorization of tensors[END_REF][START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF][START_REF] Smith | DMS: Distributed sparse tensor factorization with alternating least squares[END_REF]. The implementations in the first group are very useful tools that enable fast prototyping. Those in the second group and similar work are helpful when data fits into the memory of a single machine, which is nowadays large enough to accommodate tensors from many applications. Those in the third and fourth groups enable computations on tensors that do not fit into the memory. The ones in the third group are not designed for high performance, as MapReduce paradigm is meant to perform multiple passes over out-of-core data and perform global communication shuffling the input data.

The the best of our knowledge, there is no high-performance distributed memory implementation of algorithms for the sparse Tucker decomposition. Among the cited references above, Haten2 [START_REF] Jeon | Haten2: Billion-scale tensor decompositions[END_REF] is a MapReduce based Tucker-ALS implementation. Li et al. [START_REF] Li | An input-adaptive and in-place approach to dense tensor-times-matrix multiply[END_REF] investigate efficient shared memory execution of tensor times matrix products and as a future work mention how this can be used to perform intra-node TTM computations in a distributed memory setting.

Our data structures and the proposed reorganization of the repeated TTMc's are based on two of the cited work [6,[START_REF] Kolda | Scalable tensor decompositions for multi-aspect data mining[END_REF]. In particular, the different tensor times matrix product schemes of Kolda and Sun is used in building the TTM-tree. We make use of the data structures of Baskaran et al. [6] to represent the tensors at the nodes of the TTM-trees. Baskaran et al. also discuss a simple yet very effective scheme to increase the data reuse. This scheme computes half of the tensor matrix-chain multiplies in increasing order of modes, and the other in decreasing order of modes. While computing the first half, the factor matrices corresponding to the second half are not modified and hence one can accumulate and re-use them. A similar observation holds for the second half. Our TTM-tree representation takes this scheme one step ahead and increase data reuse at the expense of extra storage. We used PaToH [START_REF] Aykanat | PaToH: A Multilevel Hypergraph Partitioning Tool[END_REF] (version 3.2) with default options to partition the hypergraphs that we formed. We created all partitions offline, and ran our experiments on these partitioned tensors on the cluster. We do not report timings for partitioning hypergraphs with PaToH, which is quite costly, for two reasons. First, in most applications, the tensors from the real-world data are built incrementally; hence, a partition for the updated tensor can be formed by refining the partition of the previous tensor. Also, one can decompose a tensor multiple times with different ranks of approximation [START_REF] Kiers | A fast method for choosing the numbers of components in Tucker3 analysis[END_REF], potentially amortizing the time spent in partitioning across multiple runs. Second, we had to partition the tensors on a system different from the one used for the experiments. It would not be very useful to compare the runtimes from different systems and repeating the same runs on different system would not add much to the presented results. We also note that the proposed algorithms are independent from the partitioning method, and hence any hypergraph partitioning method can be used.

RR n°8801

We used the same four tensors as in previous work [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF]. These tensors are from real world data. We provide the dimensions of these tensors in Table 1. Netflix tensor has user × movie × time dimensions, which we formed from the Netflix Prize competition [START_REF] Bennett | The netflix prize[END_REF]. In this tensor, the nonzeros correspond to the user reviews of movies, and review date extends the data to the third dimension. The values of the nonzeros are determined by the corresponding review scores given by the users. We obtained The NELL-B tensor from the Never Ending Language Learning (NELL) knowledge database of the "Read the Web" project [START_REF] Carlson | Toward an architecture for never-ending language learning[END_REF], which consists of tuples of the form (entity, relation, entity) such as ('Chopin','plays musical instrument','piano'). The nonzeros of this tensor correspond to these entries discovered by NELL from the web, and the values are set to be the "belief" scores given by the learning algorithms used in NELL. Delicious and Flickr are the datasets for the web-crawl of Delicious.com and Flickr.com during 2006 and 2007, which is formed by Görlitz et al. [START_REF] Görlitz | Pints: Peer-to-peer infrastructure for tagging systems[END_REF]. These datasets consist of tuples of the form (users × resources × tags × time); hence we naturally form 4-mode tensors out of these tuples. In our experiments, we set the rank of approximation R 1 = R 2 = R 3 = 10 for 3-mode tensors, and R 1 = R 2 = R 3 = R 4 = 5 for the 4-mode tensors, ran the Tucker-ALS for 10 iteration, and reported the average time spent per Tucker-ALS iteration. We first show the effects of flat TTM-tree and the binary TTM-tree. Then we investigate the scalability.

Table 2 gives the actual runtimes (in seconds) with both flat and binary TTM-trees for comparison. The choice of TTM-tree only pertains to the local TTM operations and is independent from the distributed memory parallelism; hence we provide the runtimes for the sequential as well as parallel execution times over 2048 cores of both schemes. As seen in this table, binary TTMtree improves the runtime for all tensors during both sequential and 2048-way parallel executions. In the sequential execution, using binary TTM-tree reduced the overall execution by %27, %26, %17, and %13 for NELL-B, Netflix, Delicious, and Flickr, respectively. Similarly in the parallel execution, the overall time for using the binary TTM-tree was reduced by %8, %12, %54, and %14, respectively. These results demonstrate that even for 3 and 4-mode tensors binary TTM-tree can significantly reduce the operation count by utilizing the overlaps in the partial results of TTMS. We expect this benefit to increase even further for higher dimensional tensors.

We investigate the scalability of the proposed algorithms using flat TTM-trees in order to be able provide runtime results for runs with smaller number of processors, for which memory availability becomes a problem for the binary TTM-tree scheme. In Figure 2 we give the scalability results of our algorithm using different partitions in our dataset. For all tensors, we report the results for fine-hp and coarse-hp, which stand for using the hypergraph partitioning on fine-grain and coarse-grain task definitions, respectively. Additionally, to evaluate the hypergraph partitioning, on the Flickr tensor we provide results for random partitioning of tasks in fine-grained task scheme, which we denoted as fine-random, and block partitioning of the tasks corresponding to the consecutive rows of the U matrices, which we called coarse-block.

In our results, we first observe in Figure 2 that in all test cases our methods could scale up to 2048 cores (or 256 MPI ranks), except for fine-hp on Netflix tensor, for which we observed a slight slowdown at 2048 cores. Using the fine-hp, we were able to obtain speedups values 742x, 377x, 320x, and 135x for Flickr, Netflix, Delicious, and NELL-B tensors, respectively. This is a very promising result, particularly considering the fact that on the same dataset, in our recent work in the CP-ALS setting [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF], we could only scale to 512 cores. We observed two reasons for this relative improvement in scalability. First, our CP-ALS implementation did not employ any shared-memory parallelism, and assigned one MPI thread per core. This resulted in a high number of messages that needed to be sent/received by each MPI rank; as a result, the communication latency became a bottleneck beyond 512 cores. Here we also employ shared memory parallelism, and therefore use significantly less number of MPI-threads per node (4 vs 32), which alleviates this problem simply by employing less MPI ranks. Second, the communication pattern and the volume of CP-ALS and Tucker-ALS iterations are similar; however, Tucker-ALS involves significantly more computational load for the same amount of communication. This increase in the computational density made our algorithms less communication-bound, and more scalable. We given Table 3 to investigate computation and communication requirements with different partitions of the Flickr tensor in one Tucker-ALS iteration for 2048-way parallel run with 256 MPI ranks. In this table nzmax and nrmax refer, respectively, to the maximum tensor nonzeros per MPIrank and the maximum number of rows on which an MPI-rank has nonzeros. The first observation is that the number of tensor nonzeros per MPI-rank is always well balanced with the fine grain computations. This results in load balanced TTMc computations across MPI-ranks. On the other hand, with the coarse grain formulation, the nonzeros per MPI-rank is not well balanced. This is mostly due to feasibility; some N -1 dimensional sub-tensors (slices, when N = 3) contain much more nonzeros than others, and the coarse grain formulation cannot simply achieve better balance. The second observation is that the fine grain achieves balance on the number of rows per MPIrank as well. For example, the total number of rows in fine-hp per each mode is 130636,1450259, 28916044, and 5035767 (compare the increases with respect to the sizes in Table 1). This yields 16%, 55%, 96%, and 63% load imbalance during SVD step. Coarse-hp obtains smaller number of rows per MPI-rank (no publication) but achieves worse imbalance than the fine-hp (about 700%, 216%, 180%, 208%). As we also observed before [START_REF] Kaya | Scalable sparse tensor decompositions in distributed memory systems[END_REF], fine-hp is more effective in reducing the communication cost than the coarse-hp.

Conclusion

We discussed efficient parallelization of the alternating least squares based algorithms (Tucker-ALS also called HOOI) for Tucker decomposition of sparse tensors in current distributed memory systems. We introduced a tree-based computation scheme called TTM-tree for tensor matrix multiplication (TTM) operations which enables data reuse in the computational core of the targeted algorithms. We discussed coarse and fine grain algorithms for implementing the TTMs. The granularity in TTMs also affects the design choices for the computing the singular vectors. We discussed a careful implementation with popular libraries in which the communication cost reduction and load balance are implicitly achieved by aligning the data partition for this step with that of TTMs. We put all of these together and implemented a complete Tucker decomposition algorithm with shared (via OpenMP) and distributed memory (with MPI) parallelism. Using this implementation, we reported speedups up on up to 2048 cores.

This tensor matrix multiplication tree is of independent interest as it can also be used in sequential computing environments with dense tensors. There are many possible trees, but we explored only two of them in this study, and left further investigation as a future work. In particular, automated means to decide the shape of the tree according to a given sparse tensor need to be developed. Also, tensor times matrix-chain products are used in other algorithms [START_REF] Eldén | A Newton-Grassmann method for computing the best multilinear rank-(r 1 , r 2 , r 3 ) approximation of a tensor[END_REF] for Tucker decomposition. The proposed computational tools, including the task definitions and the TTM-tree can be used in those algorithms.

RR n°8801

Figure 2 -

 2 Figure 2 -Time for parallel Tucker-ALS iteration on Netflix and NELL-B, and speedups on Flickr and Delicious.

Table 1 -

 1 Size of tensors used in the experimentsWe conducted our experiments on the IDRIS Ada cluster, which consists of 304 nodes each having 128 GBs of memory and four 8-core Intel Sandy Bridge E5-4650 processors running at 2.7 GHz. We ran our experiments up to 2048 cores (64 nodes), which is the maximum allowed in the cluster. We tried different configurations of 32 cores within a node for the MPI/OpenMP rank/thread assignments, and found the assignment of 4 MPI ranks (and 8 OpenMP threads per MPI-rank) per node to give the best results. All codes we used in our benchmarks were compiled using the Intel C++ compiler (version 14.0.1.106) with -O3 option for compiler optimizations, -openmp flag to enable OpenMP, and -mkl option to use the Intel MKL library (version 11.1.1) for LAPACK, and BLAS routines. To obtain the sequential runtimes, we allocated 1 MPI rank with 1 threads, disabled the multithreading within MKL, and used a high-memory nodes in the cluster having 256 GB memory.

	Tensor	I 1	I 2	I 3	I 4	#nonzeros
	Netflix	480K 17K	2K	-	100M
	NELL-B 3.2M	301	638K	-	78M
	Delicious 1.4K 532K 17M 2.4M	140M
	Flickr	731	319K 28M 1.6M	112M
	5 Experimental Results

Table 2 -

 2 Comparison of the sequential and parallel run times (in secs) with fine-hp using flat and binary TTM-trees

	TTM-tree/#procs Netflix NELL-B Delicious Flickr
	flat/1	166	176	669	772
	btree/1	123	129	560	678
	flat/2048	0.44	1.3	2.09	1.04
	btree/2048	0.39	1.2	0.97	0.90

Table 3 -

 3 Statistics for the computation and communication requirements with different partitionings of the Flickr tensor in one Tucker-ALS iteration for 2048-way parallel run with 256 MPI ranks.

	Mode	Comp. load nzmax nrmax	Comm. volume Num. msg. Max Avg Max Avg
				fine-hp
	1	441K	590	2218	2029	510	509
	2	441K	8778	24K	17K	507	491
	3	441K	221K	166K	11K	265	109
	4	441K	32K	77K	53K	510	508
			fine-random
	1	443K	668	5884	2597	490	464
	2	443K	98K	409K	385K	510	510
	3	443K	545K 1744K 1735K 510	510
	4	443K	110K	432K	413K	510	510
			coarse-hp
	1	718K	22	4910	1213	255	254
	2	810K	2700	118K	66K	255	255
	3	798K	197K 3187K 810K	255	255
	4	2368K	13K	170K	102K	255	255
			coarse-block
	1	958K	252	18K	907	220	162
	2	756K	5401	126K	80K	254	252
	3	441K	130K 3324K 1250K 254	223
	4	2518K	60K	296K	138K	246	237

RR n°8801

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Inovallée 655 avenue de l'Europe Montbonnot


Clarisse Dupin

Laurent Siproudhis
email: laurent.siproudhis@chu-rennes.fr

Sébastien Henno

Sophie Minjolle

Cédric Arvieux

Pierre Tattevin

Prof Laurent

Use of human papillomavirus genotyping and biomarkers for targeted screening of anal dysplasia in human immunodeficiency virus-infected patients

Keywords: HIV, HPV, anal dysplasia, high-risk HPV genotypes, high-risk HPV E6/E7 mRNA, p16 INK4A /Ki-67

Background. Screening for anal dysplasia in human immunodeficiency virus (HIV)-infected patients is not standardized. Highresolution imaging is not adequate for mass screening, and anal cytology requires expertise. New biomarkers, selected because of their use in cervical cancer mass screening, have been originally tested for targeted and easy-to-perform screening. Methods. 120 HIV-infected individuals (males 96.4%, mean age 47±11 years) were referred for clinical examination, anoscopy, and cytological studies on anal swab. Dysplasia grading, Human Papilloma Virus genotyping, E6/E7mRNA detection and p16 INK4A /Ki-67 immunostaining were performed. High-grade lesions were histologically confirmed by anal biopsies after high-resolution anoscopy.

Results. Among the 120 anal swabs analysed, 36 (30%) had low grade, and 6 (5%) had high-grade lesions. Virus genotype was identified in 88 patients (73.3%), and 77 (64.2%) were positive for high-risk genotype(s). High-risk genotype was associated to lowgrade or high-grade lesions with a sensitivity of 0.93, and a specificity of 0.51. For E6/E7mRNA, sensitivity and specificity for lowgrade and high-grade lesions were, respectively, 0.88, and 0.78. Combination of genotyping, E6/E7mRNA and p16 INK4A /Ki-67 appropriately ruled out dysplasia in 55% of patients. Conclusions. Three routine biomarkers may avoid unnecessary invasive procedures with the perspective of an improvement of patient compliance. A decision making algorithm, based on the combination of these 3 biomarkers, is proposed.

A c c e p t e d M a n u s c r i p t 2 FUNDING Hologic (Inc. Bedford, MA, USA) offered the test kits for E6/E7 viral mRNA detection of 14 high risk HPV types (APTIMA assay). There was no other funding.

A c c e p t e d M a n u s c r i p t

INTRODUCTION

The incidence of anal cancer has dramatically increased over the last decade in high-risk groups, including men who have sex with men (MSM), women with a previous history of cervical or vulvar cancer, and immunocompromised patients, especially those who are Human Immunodeficiency Virus (HIV)-infected [START_REF] Piketty | High prevalence of anal squamous intraepithelial lesions in HIV-positive men despite the use of highly active antiretroviral therapy[END_REF][START_REF] Hoots | Human papillomavirus type distribution in anal cancer and anal intraepithelial lesions[END_REF][START_REF] Palefsky | Effect of highly active antiretroviral therapy on the natural history of anal squamous intraepithelial lesions and anal human papillomavirus infection[END_REF][START_REF] Conley | Factors associated with prevalent abnormal anal cytology in a large cohort of HIVinfected adults in the United States[END_REF]. Several series estimated that the risk of anal cancer in HIV-infected MSM is comparable to the risk of cervical cancer in the general female population [START_REF] Palefsky | Human papillomavirus-related disease in people with HIV[END_REF][START_REF] Ferlay | Estimates of cancer incidence and mortality in Europe in 2008[END_REF]. As for cervical cancers, intra-epithelial lesions induced by human papillomavirus (HPV) precede invasive anal cancer. Studies on the prevalence of HPV and its association to dysplasia in anal samples have demonstrated that both HPV and cytological abnormalities are common in high-risk populations.

The overall prevalence of HPV infection in anal carcinoma was estimated at 85-97%, HPV-16 being the most frequent genotype, found in ~75% of cases [START_REF] Vuyst | Prevalence and type distribution of human papillomavirus in carcinoma and intraepithelial neoplasia of the vulva, vagina and anus: a meta-analysis[END_REF][START_REF] Abramowitz | Human papillomavirus genotype distribution in anal cancer in France: the EDiTH V study[END_REF]. Studies based on systematic cytology of anal swabs yield a much greater burden of abnormal findings than studies based on biopsies restricted to macroscopic lesions. In a large prospective study of HIV-infected MSM, the incidence of anal intra epithelial neoplasm was 37% over a 3-year period [START_REF] De Pokomandy | HAART and progression to high-grade anal intraepithelial neoplasia in men who have sex with men and are infected with HIV[END_REF]. The prevalence of abnormal cytology ranges from 30% to 67% in these populations, and high-grade squamous intra-epithelial lesions (H-SIL) are found in 5-24%. These data advocate for the implementation of systematic screening for anal neoplasia in specific groups [START_REF] Conley | Factors associated with prevalent abnormal anal cytology in a large cohort of HIVinfected adults in the United States[END_REF][START_REF] Ciobotaru | Prevalence and risk factors for anal cytologic abnormalities and human papillomavirus infection in a rural population of HIV-infected males[END_REF]. However, screening for anal dysplasia in high-risk populations is neither universally recommended, nor standardized. French guidelines recommend annual anoscopy in HIV-infected MSM, and in HIV-infected women with HPV-related genital lesions, and encourages systematic cytological studies of anal swabs during these screening [START_REF] Morlat | Ministère des Affaires Sociales et de la Sante[END_REF], but these recommendations are poorly implemented for the following reasons: i) the number of specialized clinics offering expert evaluation and anoscopy is not sufficient in most areas; ii) this invasive procedure may be poorly tolerated by patients, and time-consuming for proctologists; iii) most pathologists are not familiar with anal cytology. In addition, both methods (i.e. anoscopy and cytological studies on anal swab) have questionable DNA genotyping and p16/Ki-67 on anal swabs had the highest sensitivity (92-100%) for the detection of anal dysplasia in high-risk patients [START_REF] Wentzensen | Human papillomavirus genotyping, human papillomavirus mRNA expression, and p16/Ki-67 cytology to detect anal cancer precursors in HIV-infected MSM[END_REF]. In order to improve the efficiency of screening programs, the aims of the present study were two-fold: i) to evaluate the diagnostic accuracy of these biomarkers, and ii) to determine their potential use for targeted screening of anal dysplasia in HIVinfected patients.

MATERIALS AND METHODS

Study design. This cross-sectional study was conducted in the Rennes University Hospital, a tertiary care teaching hospital which serves as the referral centre for HIV-infected patients, and for anal pathology in the area (Western France, catchment population estimated at one million inhabitants). All HIV-infected patients were invited for anal neoplasia screening at the department of proctology, with special emphasis in MSM and women with HPV-related genital lesions, according to national guidelines. The first visit consisted in clinical examination, digital rectal examination, anoscopy, and collection of anal swabs for cytological studies. Data about life style and HIV disease were collected using a standardized questionnaire: gender, age, sexual activity, past history of sexually transmitted infection(s), and HIV-related parameters, including CDC stage, antiretroviral treatment (ART), CD4 cell counts, and HIV viral load. Before anoscopy, two Dacron swabs were collected from the anal canal and immediately suspended into ThinPrep PreservCyt medium (Hologic, Inc. Bedford, MA, USA): one was sent for cytology and immunohistochemistry studies, the other one for virological analysis. Specimens were maintained at +4°C before processing for analyses or aliquoted and stored at -80°C. All patients diagnosed with abnormal cytology through this first step underwent high-resolution anoscopy (HRA), and guided-A c c e p t e d M a n u s c r i p t 6 biopsies were performed for histological studies if lesions were suspected, using acetic acid and Lugol's solution for better visualization. No biopsy was obtained from patient with normal cytology.

Cytology and histology. Two slides were prepared with the same medium. One was stained for cytology and the other one for immunochemistry. The pathologist who analysed samples has a long-standing experience in the analysis of anal samples, and was blinded to the results of biomarkers (HPV genotype, HPV E6/E7 mRNA and p16/Ki-67). Cytology grades were reported according to the Bethesda classification: no intraepithelial lesion (NIL), atypical squamous cell of undetermined significance (ASC-US), low-grade or high-grade squamous intraepithelial lesions (L-SIL or H-SIL) [START_REF] Solomon | The 2001 Bethesda System: terminology for reporting results of cervical cytology[END_REF]. The same pathologist classified histology as normal, low-or high-grade intra-epithelial lesions. p16 INK4A /Ki-67 dual staining. All specimens were tested for dual immunostaining p16 INK4A /Ki-67 CINtec ® PLUS (Roche, MTM laboratories AG, Heidelberg, Germany), according to the manufacturer's instructions. Protein p16 INK4A is a cell cycle regulatory protein overexpressed in high risk (hr) HPV persistent infection, and Ki-67 is a proliferation marker. The kit is a combination of two antibodies against p16 INK4A and Ki-67 used in the same slide, with different chromogens. The presence of at least one double-stained cell is interpreted as a positive test.

HPV DNA detection and genotyping. DNA testing was performed using the PreservCyt medium: 200 L were used for DNA extraction on MagNA Pure LC (Roche, Bâle, Switzerland). Anal specimens were tested for 35 HPV genotypes with CLART Human Papillomavirus 2 (Genomica, Madrid, Spain), a commercial kit for the detection of 15 low risk (lr) HPV [START_REF] Ferlay | Estimates of cancer incidence and mortality in Europe in 2008[END_REF][START_REF] Morlat | Ministère des Affaires Sociales et de la Sante[END_REF]40,42,43,44,54,61,62,71,72,81,83,84,89), and 20 high risk (hr) HPV [START_REF] Chiao | Screening HIV-infected individuals for anal cancer precursor lesions: a systematic review[END_REF][START_REF] Goldie | Cost-effectiveness of screening for anal squamous intraepithelial lesions and anal cancer in human immunodeficiency virus-negative homosexual and bisexual men[END_REF]26,31,33,35,39,45,51,52,53,56,58,59,66,68,70,73,82,85). The tests included a cellular control and an internal control. Microarray studies were analyzed using Clinical Arrays Reader High-risk (hr) HPV E6/E7 mRNA detection. Specimens were tested for E6/E7 mRNA using HPV APTIMA assay (Hologic) on the PANTHER System, which detects E6/E7 viral mRNA of 14 hr-HPV types collectively by transcription-mediated amplification [START_REF] Chiao | Screening HIV-infected individuals for anal cancer precursor lesions: a systematic review[END_REF][START_REF] Goldie | Cost-effectiveness of screening for anal squamous intraepithelial lesions and anal cancer in human immunodeficiency virus-negative homosexual and bisexual men[END_REF]31,33,35,39,45,51,52,56,58,59,66,68). An aliquot of 1 mL of each sample was transferred into the APTIMA transfer tube and introduced into the PANTHER system without any other manipulation. According to the manufacturer's instructions, a cervix sample specimen is considered positive when the value of the signal-to-cut-off ratio (S/CO) is >1. Because there is no validated S/CO threshold for anal samples, we reported the quantitative result of the hr-HPV E6/E7 mRNA test, and empirically used a conservative S/CO ratio threshold of 6. To validate this threshold, we tested HPV DNA negative samples and found that they all were either negative for hr-HPV E6/E7 mRNA, or had a S/CO ratio < 6. 1. The mean age was 47±11 years (range, 21-79).

Ethics. The database was authorized by the national regulatory institution, Commission

Among the 108 men, 89 (82.4%) were MSM. All patients underwent a clinical examination of the peri-anal skin and anal canal through an anoscope. Macroscopic signs of HPV infection were observed in 20 patients (16.7%), including external condyloma (n=12), and/or endocanal lesions (n=13). There was no ulceration and no clinical cancer. Cytological findings. Among the 120 anal swabs analyzed, 51 (42.5%) yielded NIL, 27 (22.5%) ASC-US, 36 (30%) L-SIL, and 6 (5%) H-SIL. All patients with a cytological diagnosis of H-SIL were histologically confirmed, with high-grade intra-epithelial dysplasia identified through HRA-guided biopsies. No invasive cancer was diagnosed at baseline.

HPV DNA genotyping, hr-HPV E6/E7 mRNA detection, and p16 INK4A /Ki-67 dual staining. HPV genotype was identified in 88 patients (74.6%), no HPV was found in 30 patients (25.4%), and HPV was detected, but could not be genotyped in the remaining two patients (both had NIL). Table 2 highlights the prevalence of high risk genotypes according to the cytological status. The following genotypes were the most common that differed according to cytological status: 16, 31, 35, 51, 52. Table 3 summarizes the mean number of different genotypes and the prevalence of hr-HPV DNA according to cytological status. The number of hr-HPV genotypes was significantly associated with Bethesda grading. Among specimens positive for at least one hr-HPV type, the prevalence of multiple infections was 75.3%, with a mean number of 3.5 ± 2.4 different genotypes (range, 2-11). Hr-HPV DNA was significantly associated with cytological grade, as were positivity of hr-HPV E6/E7 mRNA, mean S/CO value, and p16 INK4A /Ki-67 dual staining. Due to sub-optimal quality of samples slides and lack of material, this latter test could be performed in only 74 patients, and was not interpretable in 10 of these. The diagnosis accuracy of the three tests has been highlighted in the Table 4. For the diagnosis of L-SIL or H-SIL, the identification of at least one hr-DNA genotype exhibited a sensitivity (Se) of 0.93, a specificity (Sp) of 0.51, a negative predictive value (NPV) of 0.88, and a positive predictive value (PPV) of 0.66. High risk HPV E6/E7 mRNA assay had Se, Sp, NPV, and PPV of, respectively, 0.88, 0.78, 0.87, and 0.80. For p16 INK4A /Ki-67 testing, Se, Sp, NPV, and PPV were 0.64, 0.90, 0.68, and 0.88, respectively. Among the 41 patients with no hr-HPV (34.7%), 3 (7.3%) had L-SIL, but none had H-SIL. Of the 23 patients (31.1%) with positive p16 INK4A /Ki-67 testing, 16 (69.6%) had L-SIL, and 3 (13%) had H-SIL. Fig. 1 depicts the triage of our population based on the results of these tests. Briefly, a negative result on hr-HPV E6/E7 mRNA test could have avoided subsequent evaluation in one half of the population study, as 91.2% of them had no dysplasia and none had H-SIL. On the other hand, in patients positive for hr-HPV E6/E7 mRNA and p16 INK4A /Ki-67 tests (18.6% of the cohort), a careful examination is required as >85% of them presented with L-SIL or H-SIL. The same would apply to the 14.4% of patients with negative p16 INK4A /Ki-67 but carrying HPV-16 or HPV-18 DNA, as this combination was associated to dysplasia in ~50% of cases. A last population at low risk of anal dysplasia was identified in this cohort: patients with negative p16 INK4A /Ki-67, no HPV-16 or HPV-18, and <3 hr-HPV genotypes: only one out of 8 had dysplasia (L-SIL). To summarize, careful interpretation of hr-HPV E6/E7 mRNA, hr-HPV genotyping and p16 INK4A /Ki-67 tests in this cohort could accurately classify the risk of dysplasia in 90.8% of patients (Fig. 1).

DISCUSSION

The increasing incidence of anal cancer in HIV-infected patients, often diagnosed at a late stage, advocates for the implementation of systematic or targeted screening of high-risk patients. However, no strategy for anal cancer screening has been validated thus A c c e p t e d M a n u s c r i p t 10 far in this population. Hence, practices vary widely from one cohort to another, as a function of patients will, scientific interest of the physicians in charge, and availability of specialized centres qualified for anal cancer screening. In addition, although HRA with macroscopically-guided biopsies is highly specific, its sensitivity was recently challenged in a study comparing HRA with a combination of different cytological testing [START_REF] Wentzensen | Human papillomavirus genotyping, human papillomavirus mRNA expression, and p16/Ki-67 cytology to detect anal cancer precursors in HIV-infected MSM[END_REF]. Lastly, anal cytology requires pathologist expertise: Inter-agreement assessment has been recently found to be moderate-to-good in HIV-infected MSM [START_REF] Darragh | Interrater agreement of anal cytology[END_REF], and the sensitivity and specificity of anal cytology has been shown to be 70-93% and 67-92%, respectively, as compared with histological studies from biopsies [START_REF] Fox | The value of anal cytology and human papillomavirus typing in the detection of anal intraepithelial neoplasia: a review of cases from an anoscopy clinic[END_REF][START_REF] Chiao | Screening HIV-infected individuals for anal cancer precursor lesions: a systematic review[END_REF][START_REF] Nathan | Performance of anal cytology in a clinical setting when measured against histology and high-resolution anoscopy findings[END_REF]. For H-SIL, NPV of anal cytology were excellent (mean, 85%; range 76-92) in Nathan et al. series. Of note, these authors have shown that the sensitivity of anal swab increases with the extension of the disease, and reaches 86% when at least 2 quadrants are involved [START_REF] Nathan | Performance of anal cytology in a clinical setting when measured against histology and high-resolution anoscopy findings[END_REF].

Referral to specialized centres performing HRA and cytology may be cost-effective in high-risk population, especially HIVinfected MSM and HIV-infected women with HPV-related genital lesions [START_REF] Goldie | Cost-effectiveness of screening for anal squamous intraepithelial lesions and anal cancer in human immunodeficiency virus-negative homosexual and bisexual men[END_REF]. However, reliable tests that would allow a preselection of this high-risk population would be most welcome, by decreasing the financial costs of systematic screening, the workload for specialized centres, and the number of patients undergoing these invasive, uncomfortable tests. On the one hand, biomarkers are costly (approximately 35 to 50 Euros per test); on the other, HRA is somewhat a long procedure (30-45 minutes), the assessment of elementary lesions needs a learning curve, and the material is costly (6000 to 18000 Euros).

Our study suggests that the combination of three biomarkers, hr-HPV E6/E7 mRNA, hr-HPV genotypes and p16 INK4A /Ki-67, would allow targeted screening of anal dysplasia in HIV-infected patients, by ruling out L-SIL and H-SIL in > 55% of patients. It must be outlined that no case of H-SIL would have been overlooked with a strategy that would perform anal cytology and histology only in patients positive with either hr-HPV E6/E7 mRNA and/or p16 INK4A /Ki-67. By contrast, the guided-biopsies were performed for histological studies if lesions were suspected, using acetic acid and Lugol's solution for better visualization. According to ethical We speculated that defining viral profiles and sensitive biomarkers might help to develop diagnostic strategies and tailored screening in high-risk groups. We selected candidate biomarkers of viral oncogenicity (hr-HPV DNA, hr-HPV E6/E7 mRNA) and cytological abnormalities (p16 INK4A /Ki-67), which has already proved useful in screening of cervical cancer. They all can be routinely performed and automatized in a multicentre setting. More than 100 HPV genotypes have been linked to anal carcinoma, with several high-risk subtypes identified, HPV-16 being the most common, although the proportion of non-HPV-16 oncogenic genotypes has increased in recent series [START_REF] Wong | Human papillomavirus genotypes in anal intraepithelial neoplasia and anal carcinoma as detected in tissue biopsies[END_REF]. In our study, HPV-16 was found in 35% of swabs positive for hr-HPV, and 6 other oncogenic genotypes were found in > 10% of patients: HPV 52, 58, 66, 33, 53, and 51. It emphasizes the need of genomic tests which detects all hr-HPV genotypes: the APTIMA assay is one of them. For screening purposes, biomarkers must primarily be sensitive, if they are to be used as 'rule-out' tests. In this regard, the detection of hr-HPV E6/E7 mRNA appears as a suitable test: in the study presented herein, one-half of the population had no detection of hr-HPV E6/E7 mRNA, and 91.2% of them had no dysplasia (neither L-SIL nor H-SIL). On the contrary, p16 INK4A /Ki-67 would not qualify as a good biomarker to rule out anal dysplasia, given its low sensitivity (0.64), but would be of interest as a confirmatory test, given its high specificity: only 7.7% of patients with positive p16 INK4A /Ki-67 had normal cytology. However, not all patients were tested because the material was no longer available when p16INK4A/Ki67 was retrospectively performed. These preliminary results must be confirmed by a prospective study.

Using three reproducible tests (hr-HPV E6/E7 mRNA, hr-HPV genotype and p16 INK4A /Ki-67) with no need for expertise, we selected 53 patients (44.9%) who would be candidates for pathologist examination of anal swabs (Fig. 1). Of these, 36 (67.9%) had L-SIL or H-SIL on anal cytology and 9 (17%) presented dysplasia on subsequent examination, indicating the need for careful A c c e p t e d M a n u s c r i p t 12 examination of the anal canal with guided biopsies. For the other half of this high-risk HIV-infected population, hr-HPV E6/E7 mRNA, and hr-HPV genotype could rule out anal dysplasia, and avoid further tests. The diagnostic test performance using a combined analysis was lacking in the present study. The aim of the present study was in fact to evaluate the accuracy of each biomarker, and the size of the study population does not allow the calculation of predictive values among small subgroups. We tried to design a simple algorithm in figure 2 which remains to be tested prospectively among a larger cohort in order to analyse the performance of a stepwise combination of tests.

This study has limitations, especially in the lack of systematic follow-up data in patients with ASC-US or L-SIL. However, studies of HPV clearance, co-infections or replacement by new HPV DNA subtypes are in progress in our centre. Changes in cytological status over time may be the result of subsequent infectious events, thus the evolution of cytological patterns remains speculative without concomitant follow-up of HPV status. Monitoring the rate of progression of various stages of anal dysplasia in a large HIV-positive cohort, de Pokomandy et al. found no regression of extensive untreated dysplasia, and a dramatic increase in the cumulative proportion of subjects with L-SIL or H-SIL at any time point since enrolment in the cohort [START_REF] De Pokomandy | HAART and progression to high-grade anal intraepithelial neoplasia in men who have sex with men and are infected with HIV[END_REF]. These data suggest that HPV clearance is not a frequent event in HIV-infected patients. Other limitations of our study include its monocentric design, and the absence of a validation cohort, as diagnostic algorithms tend to be less accurate when applied in other cohorts. In addition, as HRA and biopsies were only performed in patients with dysplasia, histological studies were not performed in the whole cohort.

However, to the best of our knowledge, only few studies on the yield of biomarkers for the screening of anal canal cancer in high-risk groups have been published [START_REF] Wentzensen | Human papillomavirus genotyping, human papillomavirus mRNA expression, and p16/Ki-67 cytology to detect anal cancer precursors in HIV-infected MSM[END_REF][START_REF] Salit | Screening for HIV-associated anal cancer: correlation of HPV genotypes, p16, and E6 transcripts with anal pathology[END_REF][START_REF] Silling | Human papillomavirus oncogene mRNA testing for the detection of anal dysplasia in HIV-positive men who have sex with men[END_REF][START_REF] Phanuphak | Use of human papillomavirus DNA, E6/E7 mRNA, and p16 immunocytochemistry to detect and predict anal high-grade squamous intraepithelial lesions in HIV-positive and HIVnegative men who have sex with men[END_REF], and only one studied high-risk HPV E6/E7 mRNA and p16 INK4A /Ki-67 immunostaining simultaneously [START_REF] Wentzensen | Human papillomavirus genotyping, human papillomavirus mRNA expression, and p16/Ki-67 cytology to detect anal cancer precursors in HIV-infected MSM[END_REF]. All four studies published to date confirm that these biomarkers may be of interest for anal dysplasia screening, although further investigations are required to better define their respective added values. The study presented A c c e p t e d M a n u s c r i p t 13 herein is the second where hr-HPV E6/E7 mRNA was detected with APTIMA assay, which detects 14 hr-HPV, and has been associated with improved diagnostic yield for cervical dysplasia [START_REF] Szarewski | Comparison of predictors for high-grade cervical intraepithelial neoplasia in women with abnormal smears[END_REF][START_REF] Burger | HPV mRNA tests for the detection of cervical intraepithelial neoplasia: a systematic review[END_REF]. In summary, easy to perform, reproducible and non-invasive biomarkers, such as the identification of high-risk HPV E6/E7 mRNA, hr-HPV genotypes, and p16 INK4A /Ki-67 immunostaining on anal swab may be used for the screening of anal dysplasia in HIV-infected patients. In our study, their combination appropriately ruled out L-SIL and H-SIL in >55% of patients. A simple algorithm of a stepwise testing (figure 2) remains to be tested prospectively among a larger cohort. The cost effectiveness strategies and the patient compliance need to be clarified in population studies.

  A c c e p t e d M a n u s c r i p t 5 reproducibility and sensitivity to detect dysplasia of the anal canal. Testing several biomarkers, Wentzensen et al. found that HPV

A c c e p t e d M a n u s c r i p t 7 (

 7 Genomica). A sample of 2 mL of each liquid-based medium was then frozen at -80°C. A hr-HPV DNA test was considered positive if at least one hr-HPV genotype was detected.

  Nationale Informatique et Liberté (CNIL n°1412467), and the study was approved by our institutional review board (N09-14). All participants provided written informed consent.Statistical analysis.Mean and standard deviation were used to describe quantitative variables; number and percentage were used to describe qualitative variables. Differences between groups were measured using non-parametric tests where appropriate (Wilcoxon, Pearson, Chi2 or MacNemar). A P value <0.05 was considered significant. The performance of p16 INK4A /Ki-67 dual staining, hr-HPV DNA, and hr-HPV E6/E7 mRNA for the diagnosis of L-SIL and H-SIL was estimated. Statistical analyses were performed using JMP ® Pro 9.0.2, 2010 (SAS Institute Inc., Cary, North Carolina, USA).

  A c c e p t e d M a n u s c r i p t 8 RESULTS Study participants. Between October 2010 and November 2011, 120 HIV-infected patients were enrolled in this cross-sectional study. The main characteristics of the study population are summarized in Table

  A c c e p t e d M a n u s c r i p t 11 concerns, no biopsy was obtained from patients with normal cytology: in fact, biopsies are somewhat painful and they carry a bleeding risk.

18 FIG 1 :

 181 FIG 1: Triage of human immunodeficiency virus-infected patients screened for anal dysplasia using a combination of 3

u s c r i p t 22 Hr

 22 HPV High-risk Human Papilloma Virus genotypes are classified according to their prevalence in the study. * Cytology patterns significantly different with p<0.05; NIL, no intraepithelial lesion; ASC-US, atypical squamous cell of undetermined significance; L-SIL, low-grade squamous intraepithelial lesions; H-SIL, high-grade squamous intraepithelial lesion A c c e p t e d M a n u s c r i p t 23

A c c e p t e d M a n u s c r i p

Figure 2

 2 Figure 2

TABLE 1

 1 Characteristics of enrolled patients (N=120)

	Characteristics	Number of patients (%)
	Male Gender	108 (90)
	Age (years)
	20-29	6 (5)
	30-39	19 (16)
	40-49	51 (43)
	50-59	26 (21)
	>60	18 (15)
	HIV infection stage
	A	72 (60)
	B	19 (15.8)
	C	28 (23.4)
	Undetermined	1 (0.8)
	Current cART
	Yes	112 (93.3)
	No	5 (4.2)
	Unknown	3 (2.5)
	Nadir of CD4 cell count (cell/mm3)
	<200	55 (45.8)
	200-400	42 (35)
	>400	22 (18.4)
	Undetermined	1 (0.8)
	HIV load at enrolment (copies/ml)
	<40	93 (77.5)
	40-100	11 (9.2)
	>100	15 (12.5)
	Undetermined	1 (0.8)
	Smoking
	Current and/or past	70 (58.3)
	Never	48 (40)

TABLE 2

 2 High-risk Human Papilloma Virus genotypes according to cytological status, number (%)

	Hr-HPV	NIL	ASC-US	L-SIL	H-SIL	Total
	genotype	(%, n=49)	(%, n=27)	(%, n=36)	(%, n=6)	(%, n=118)
	16*	8(16)	2(7.4)	16(44.4)	5(83.3)	31(25.8)
	52*	3(6)	4(14.8)	8(22.2)	3(50)	18(15)
	58	6(12)	1(3.7)	10(27.7)	1(16.6)	18(15)
	66	4(8)	7(25.9)	5(13.8)	1(16.6)	17(14.2)
	33	4(8)	2(7.4)	6(16.6)	2(33.3)	14(11.7)
	53	4(8)	1(3.7)	8(22.2)	1(16.6)	14(11.7)
	51*	1(2)	2(7.4)	10(27.7)	0(0)	13(10.8)
	70	3(6)	3(11.1)	4(11.1)	2(33.3)	12(10)
	35*	1(2)	2(7.4)	6(16.6)	2(33.3)	11(9.2)
	31*	1(2)	1(3.7)	7(19.4)	1(16.6)	10(8.3)
	59	3(6)	0(0)	6(16.6)	1(16.6)	10(8.3)
	18	3(6)	3(11.1)	2(5.5)	0(0)	8(6.7)
	56	3(6)	1(3.7)	2(5.5)	1(16.6)	7(5.8)
	39	3(6)	1(3.7)	2(5.5)	0(0)	6(5)
	68	1(2)	2(7.4)	2(5.5)	1(16.6)	6(5)
	82	0(0)	2(7.4)	2(5.5)	1(16.6)	5(4.

TABLE 3

 3 High-risk Human Papilloma Virus genotypes and biomarkers according to cytological status

	NIL	ASC-US	L-SIL	H-SIL	Total
	n=49	n=27	n=36	n=6	n=118

Cytology patterns are significantly different for each parameter (P<0.0001) hr

-HPV, High-risk Human Papilloma Virus; NIL, no intraepithelial lesion; ASC-US, atypical squamous cell of undetermined significance; L-SIL, low-grade squamous intraepithelial lesions; H-SIL, high-grade squamous intraepithelial lesion; S/CO, signal-to-cut-off ratio

TABLE 4

 4 Diagnostic accuracy of biomarkers in the diagnosis of anal dysplasia

	Sensitivity	Specificity	Negative	Positive
			Predictive	Predictive
			Value	Value

ACKNOWLEDGEMENTS

We are indebted to all the patients who participated in the study, and the health care workers who took care of them.