In most circumstances, you should try to avoid computing the Gram matrix A^T A at all, as it squares the condition number of A. See also my explanation in another thread on least-squares problems.
I need a sparse factorization of A^T A which would be used in solving a system of the form (A^T A)^{-1}b for many possible vectors b.
Are you solving least-square problems? That is, are you actually solving (A^T A)^{-1} A^T b? In that case what you want is to use the A=QR factorization and solve R^{-1} (Q^* b), which is equivalent but avoids squaring the condition number (or multiplying large matrices together). In fact, the QR factorization object in Julia will do this for you:
QR = qr(A) # QR factorization (sparse if A is sparse)
x = QR \ b # least-square solution
Is there some other reason why you want to factorize the Gram matrix?