Compressed image sensing in DNA sequencing

UNFINISHED WORK published only for collaboration with others

I was recently introduced to compressed sensing, so I am going to teach myself how to use it to my problems.  Instead of struggling with CS math treaties, I will go through as many programming examples as possible, like this one aired on NPR, which makes a point that piano music can be sparsely coded and decoded with a basis consisting of piano notes, but not (very well) human speech.  I liked this tutorial because it taught me about taking advantage of duality when designing the sampling strategy: if the signal is sparse in time, sample randomly in frequency; conversely, if the signal is sparse in frequency space, sample randomly in time.

Introduction

Imaging in next gen sequencing

The next-gen DNA sequencing derives its high thoughput by observing as many randomly fragmented DNA strands in a given reaction area as possible. But since DNA strand is a molecule (polymer, more precisely), it is damned near impossible to see (which is where the 3rd gen sequencing comes in), a DNA strand is amplified to a cluster big enough to be seen with high resolution collection optics, as shown in the cartoon schematic of a Solexa/Illumina cluster:

If all goes well, a high resolution camera can observe multiple clusters fluorescing with a given dye.  By aligning the images for different dyes, a base calling algorithm can classify a given cluster to be fluorescing with one of the possible dyes at a given cycle, as shown in this example, where we can manually call the bases for 2 clusters (circled):

Of course, the reality is never this clean.  To observe as many clusters in 1 frame as possible, we want to size the cluster to span as few pixels as we can get away with, as in this example of a portion of "T dye" image of an early Solexa run rendered in ImageJ:

An average cluster's PSF in these runs is on the order of a micron, so a pixel is on the order of a 200 nm.

Problem: high resolution imaging is expensive and difficult, but low resolution imaging cannot do the job

Focusing and aligning optics system down to an order of 100 nm requires many things to work exactly as they should.  Suppose that I have an imaging system that can resolve 1 um (which is not easy either!).  Even if I could focus the clusters to this poorer resolution camera (I can't), the clusters would be so smeared (binned) that I cannot pick up a putative cluster from the resulting image like this:

I want to see if the original image is sparse enough to benefit from compressed sensing.

Simulating random imaging in Matlab

Single pixel camera experiment showed that a sparse (in space) image can be reconstructed from randomly oversampled single-photo-detector data.  The random sampling was done with a DLP, which had sufficient spatial resolution.  In Matlab, taking ONE random sample means to convolve the signals with a randomly chosen basis vector--and DFT (discrete Fourier transform) is as good a basis as any.

1D DFT

To make sure I understand the operation, let's fall back to the definition of a DFT given a vector x[k], k = 0 through N-1:

X[n] = SUM_k { exp(2i n k/N) x[k] },  k and n = 0~N-1

     = [exp(2i n 0/N) ... exp(2i n (N-1)/N)] * [x[0] ... x[N-1]]'

So X can be computed with a matrix and vector multiplication: X = A x

[ X[0]  ]   [...       ...        exp(2i 0 (N-1)/N)]   [ x[0]   ]

[ X[n]  ] = [... exp(2i n k/N)   exp(2i n (N-1)/N)] * [ ...    ]

[ X[N-1]]   [... exp(2i(N-1)k/N) exp(2i (N-1)^2/N)]   [ x[N-1] ]

X = dftmtx(N) * x

The square matrix A matrix above is the DFT matrix, obtained in Matlab with the function dftmtx(N).  To invert the transform, complex conjugate of the matrix is required, because the definition of the inverse transform is identical to the forward transform, except for the negative sign in the exp argument.

x = conj(dftmtx(N)) * X .* (1/N)

An image is commonly represented as a 2D matrix, but is in fact commonly stored as 1D vector in memory.  So 1D can be used without loss of generality.  BUT, it is difficult for people to "see" an image unrolled to 1D, so a 2D treatment is useful.

Random sampling in 1D frequency space

In the Matlab CS tutorial mentioned above, random sampling in discrete frequency space is accomplished by randomly selecting K rows of the DFT matrix (which is square), which is extremely easy to do in Matlab:

qN = randperm(N); randN = dftN(qN(1:K), :);

y = randN * truth + v;

where "truth" is the unobservable true state that we are trying to estimate, and v is an additive noise term that is bound to be present all real world measurements--usually set to 0 at first pass through the theory.  Actually implementing such sensing hardware is extremely difficult: not the noise injection part--which the nature does an excellent job of--but the construction of a hardware that will pseudo-randomly mix the spectral components of the true states into under-sampled measurements.

2D DFT

Falling back on the definition of 2D DFT:

X[n,m] = SUM_k SUM_l { exp(2i n k/N) exp(2i m l/M) x[k,l] },

         where k and n = 0~N-1, l and m = 0~M-1

       = SUM_k { exp(2i n k/N) SUM_l { exp(2i m l/M) x[k,l] } }

Using a similar manipulation as the 1D case, X, which is now an NxM matrix, can be obtained like this:

[ X[  0,:] ]   [...       ...        ... exp(2i 0 (N-1)/N)]

[ X[  n,:] ] = [... exp(2i n k/N)   ... exp(2i n (N-1)/N)]

[ X[N-1,:] ]   [... exp(2i(N-1)k/N) ... exp(2i (N-1)^2/N)]

               [x[  0,0]      ...       x[  0,M-1]] 

             * [x[  n,0] ... x[n,m] ... x[  n,M-1]]

               [x[N-1,0]      ...       x[N-1,M-1]]

               [...       ...          ... exp(2i 0 (M-1)/M)]

             * [... exp(2i m l/M)     ... exp(2i m (M-1)/M)]

               [... exp(2i m (M-1)/M) ... exp(2i (M-1)^2/M)]

More concisely,

X = dftmtx(N) * x * dftmtx(M) = ((x * dftmtx(M))' * dftmtx(N)')' 

x = conj(dftmtx(N)) * X * conj(dftmtx(M)) ./ (N*M)

It seems that the concise notation can be generalized to a recursive relationship between adjacent dimension:

X[i] = (X[i-1]' * dftmtx(Ni)')', X[0] = x

X and x are multi-dimensional matrices with cardinality (N1, ... Nd).

As a sanity check, I should be able to take a round-trip to the frequency space and back:

img_orig = imread('/mnt/work/Sensing/L004/C1.1/s_4_61_t.tif');

img_orig = img_orig(1:1024, 1:1024);

img_data = double(img_orig);

imshow(img_orig, []); %[] is necessary to actually show the image

fft_orig = fft2(img_orig); imshow(real(fft_orig), []);

img_back = uint16(ifft2(fft_orig));

[N, M] = size(img_orig)

dftN = dftmtx(N); dftM = dftmtx(M);

dft_orig = dftmtx(N) * img_data * dftmtx(M); imshow(real(dft_orig), [])

img_back = uint16(conj(dftN) * dft_orig * conj(dftM) .* (1/(N*M)));

max(max(img_orig - img_back))

The result of the round trip is the maximum value of the difference between the original image and the round trip image, which turns out to be 0.

Flattening an image to 1D vector

For compressive sensing, 2D data is unwieldy, because the literature and examples all use 1D data.  A silver lining is that in computer memory, multi-dimensional memory is already flattened out to a vector, and flattening a matrix to a vector in Matlab is trivial.

img_flat = img_orig(:); % Flatten an image into a vector

 BUT flattening a 2D DFT to a 1D will NOT yield the same spectra as the 1D transform of the flattened image:

fft2D = fft2(img_orig); fft2D_flat = fft2D(:);
fft1 = fft(img_flat);
>> max(fft2D_flat - fft1)
ans =
  -8.7022e+07 + 2.7322e+07i

Sparse random sampling of image

The example image I am using is 1024x1024 pixels, with each pixel representing roughly 200 nm squared surface area.  Suppose I get a sensor with 1.6 um resolution (8 times worse than the original).

WIP