Closed liamgd closed 1 year ago
Various implementations of region mapping should be tested for performance.
Creation of region filtration boolean matrix
Filtering tensor using boolean matrix
Test of speed of creation of region filtration matrix with 300 samples, 50,000 methylation loci, and 10,000 gene expression loci.
Line # Hits Time Per Hit % Time Line Contents
...
77 # Chromosomal filtration
78 1 153939.6 153939.6 1.4 M_chrom_t[:, None] == G_chrom_t
79 1 136644.3 136644.3 1.3 G_chrom_t[:, None] == M_chrom_t
80
81 # Position filtration
82 1 2145274.1 2145274.1 19.8 torch.isclose(M_pos_t[:, None], G_pos_t, rtol=0, atol=window)
83 1 2169992.9 2169992.9 20.0 torch.isclose(G_pos_t[:, None], M_pos_t, rtol=0, atol=window)
84 1 720406.6 720406.6 6.6 (M_pos_t[:, None] < G_pos_t + window).logical_and(M_pos_t[:, None] > G_pos_t - window)
85 1 671274.2 671274.2 6.2 (G_pos_t[:, None] < M_pos_t + window).logical_and(G_pos_t[:, None] > M_pos_t - window)
86 1 1111070.2 1111070.2 10.3 (M_pos_t[:, None] - G_pos_t).abs() < window
87 1 1208527.4 1208527.4 11.2 (G_pos_t[:, None] - M_pos_t).abs() < window
88 1 1144775.9 1144775.9 10.6 ((t := (M_pos_t[:, None] - G_pos_t)) < window).logical_and(t > -window)
89 1 1358843.5 1358843.5 12.5 ((t := (G_pos_t[:, None] - M_pos_t)) < window).logical_and(t > -window)
From this, we gather that the fastest chromosomal filtration is:
G_chrom_t[:, None] == M_chrom_t
And the fastest chromStart position filtration within a window is:
(G_pos_t[:, None] < M_pos_t + window).logical_and(G_pos_t[:, None] > M_pos_t - window)
Tests of speed of combining these three boolean matrices with with 300 samples, 50,000 methylation loci, and 10,000 gene expression loci.
81 1 812456.0 812456.0 30.3 (G_chrom_t[:, None] == M_chrom_t).logical_and(G_pos_t[:, None] < M_pos_t + window).logical_and(G_pos_t[:, None] > M_pos_t - window)
82 1 895275.5 895275.5 33.4 (G_chrom_t[:, None] == M_chrom_t) + (G_pos_t[:, None] < M_pos_t + window) + (G_pos_t[:, None] > M_pos_t - window)
83 1 957611.4 957611.4 35.8 torch.stack((G_chrom_t[:, None] == M_chrom_t, G_pos_t[:, None] < M_pos_t + window, G_pos_t[:, None] > M_pos_t - window)).sum(dim=0, dtype=bool)
A series of logical_and
s is the fastest:
(G_chrom_t[:, None] == M_chrom_t).logical_and(G_pos_t[:, None] < M_pos_t + window).logical_and(G_pos_t[:, None] > M_pos_t - window)
(This example snippet is for a cis analysis, but is very similar for distal. Trans analyses would include only the first expression before the first logical_and)
Current region filtration:
if p_thresh is not None:
if region == 'all':
indices = P[:, 1] <= p_thresh
elif region == 'cis':
indices = (
(P[:, 1] <= p_thresh)
.logical_and(M_chrom_t[index - 1, None] == G_chrom_t)
.logical_and(
M_pos_t[index - 1, None] < G_pos_t + window
)
.logical_and(
M_pos_t[index - 1, None] > G_pos_t - window
)
)
elif region == 'distal':
indices = (
(P[:, 1] <= p_thresh)
.logical_and(M_chrom_t[index - 1, None] == G_chrom_t)
.logical_and(
(
M_pos_t[index - 1, None] < G_pos_t - window
).logical_or(
M_pos_t[index - 1, None] > G_pos_t + window
)
)
)
elif region == 'trans':
indices = (P[:, 1] <= p_thresh).logical_and(
M_chrom_t[index - 1, None] != G_chrom_t
)
else:
if region == 'cis':
indices = (
(M_chrom_t[index - 1, None] == G_chrom_t)
.logical_and(
M_pos_t[index - 1, None] < G_pos_t + window
)
.logical_and(
M_pos_t[index - 1, None] > G_pos_t - window
)
)
elif region == 'distal':
indices = (
M_chrom_t[index - 1, None] == G_chrom_t
).logical_and(
(
M_pos_t[index - 1, None] < G_pos_t - window
).logical_or(
M_pos_t[index - 1, None] > G_pos_t + window
)
)
elif region == 'trans':
indices = M_chrom_t[index - 1, None] != G_chrom_t
if indices_list is not None:
indices_list.append(indices)
P = P[indices]
output_sizes.append(len(P))
results.append(P)
For example, the code would be shorter if filtration tensors were generated for p-value filtration, chromosome filtration, and location window filtration. Then, element-wise boolean algebra dependent on the region and p-value threshold would determine the indices. Unfortunately, this seems to be much slower than simply enumerating every combination of region and p-value.
Times for tests with 300 samples, 50000 methylation loci, and 10000 gene expression loci on a 2070 super with 8 GB VRAM (constant calculation and loop time, saving included only for chunked calculations):
No chunks, no saving:
tecpg run mlr-full: *Too much GPU memory allocation*
tecpg run mlr-full --cis: 93.1929 seconds
tecpg run mlr-full --distal: 94.3363 seconds
tecpg run mlr-full --trans: *Too much GPU memory allocation*
tecpg run mlr-full -p 0.05: 93.2036 seconds
tecpg run mlr-full -p 0.05 --cis: 93.427 seconds
tecpg run mlr-full -p 0.05 --distal: 93.5236 seconds
tecpg run mlr-full -p 0.05 --trans: 93.5954 seconds
Chunks of 2000 methylation loci, saving included:
tecpg run mlr-full -l 2000: 169.2086 seconds
tecpg run mlr-full -l 2000 --cis: 102.665 seconds
tecpg run mlr-full -l 2000 --distal: 104.6743 seconds
tecpg run mlr-full -l 2000 --trans: 180.5167 seconds
tecpg run mlr-full -l 2000 -p 0.05: 101.5071 seconds
tecpg run mlr-full -l 2000 -p 0.05 --cis: 97.5001 seconds
tecpg run mlr-full -l 2000 -p 0.05 --distal: 98.7852 seconds
tecpg run mlr-full -l 2000 -p 0.05 --trans: 101.6967 seconds
Judging by these speeds, it seems that the current implementation of region filtration slightly lengthens run time but slightly reduces the time it takes to save (judging by the chunked analyses that included saving). It may be possible to reduce the speed further for cis, distal, and trans filtration by filtering the necessary tensors before each computation on a particular methylation locus.
This would entail creating a region filtration mask tensor before computation. The Y (or M_row if it is faster), XtXi_Xt, and X tensors would be filtered with this mask. Then, the p-value filtration mask would be calculated after and would further filter the region mask.
With the above performance improvements, we obtain the following performance measurements.
Comparison to the previous region filtration implementation to the new one. The first value in seconds is the previous time, and the second value is the new time.
No chunks, no saving:
tecpg run mlr-full: *Too much GPU memory allocation*
tecpg run mlr-full --cis: 93.1929 seconds to 41.4291 seconds
tecpg run mlr-full --distal: 94.3363 seconds to 43.6176 seconds
tecpg run mlr-full --trans: *Too much GPU memory allocation*
tecpg run mlr-full -p 0.05: 93.2036 seconds to 96.6992 seconds
tecpg run mlr-full -p 0.05 --cis: 93.427 seconds to 50.5002 seconds
tecpg run mlr-full -p 0.05 --distal: 93.5236 seconds to 45.169 seconds
tecpg run mlr-full -p 0.05 --trans: 93.5954 seconds to 133.2633 seconds
Chunks of 2000 methylation loci, saving included
tecpg run mlr-full -l 2000: 169.2086 seconds to 173.3093 seconds
tecpg run mlr-full -l 2000 --cis: 102.665 seconds to 49.2062 seconds
tecpg run mlr-full -l 2000 --distal: 104.6743 seconds to 59.2901 seconds
tecpg run mlr-full -l 2000 --trans: 180.5167 seconds to 235.1578 seconds
tecpg run mlr-full -l 2000 -p 0.05: 101.5071 seconds to 101.6262 seconds
tecpg run mlr-full -l 2000 -p 0.05 --cis: 97.5001 seconds to 47.2857 seconds
tecpg run mlr-full -l 2000 -p 0.05 --distal: 98.7852 seconds to 50.1013 seconds
tecpg run mlr-full -l 2000 -p 0.05 --trans: 101.6967 seconds to 145.8139 seconds
In my opinion, for large input matrices, the new way is faster. I believe this algorithm scales better with larger inputs because with respect to the gene expression count, the time complexity with the previous algorithm is superlinear for the computation section and the region mask creation and filtration runs in linear time. With the new algorithm, the linear mask creation and filtration happen before the superlinear computation. Assuming g gene expression loci, the number of gene expression loci for the computation would be a factor between 0 and 1 multiplied by that number. The superlinear function in which the computation runs would then scale that factor even lower, which accounts for the theoretical speed improvement of the new function for large inputs.
Moved filtration before computation in e8fc49a.
Further optimization warrants a new issue.
The regression full algorithm should have region mapping capabilities as outlined in #7.
Optional annotation files M_annot and G_annot will supplement the M and G dataframes and determine the distance between loci. They will be required for cis, distal, and trans analyses.
All: regress every methylation and gene expression pair. Cis: regress every methylation and gene expression pair within a default window on the same chromosome. Distal: regress every methylation and gene expression pair outside of a default window on the same chromosome. Trans: regress every methylation and gene expression pair on different chromosomes.