diff --git a/docs/content/counts.md b/docs/content/counts.md index 1311db60..546a430b 100644 --- a/docs/content/counts.md +++ b/docs/content/counts.md @@ -7,46 +7,45 @@ $ dnmtools counts [OPTIONS] -c ## Description -The `counts` command takes the mapped reads and produces the -methylation level at each genomic cytosine (both strands), with the -option to produce only levels for CpG-context cytosines. While most -DNA methylation in vertebrates is in the CpG context, cytosines in -other sequence contexts, such as CXG or CHH (where H denotes adenines, -thymines, or cytosines and X denotes adenines or thymines) may also be -methylated. Non-CpG methylation occurs frequently in plant genomes and -pluripotent mammalian cells such as embryonic stem cells. And possibly -whatever cells you are studying. The output of `counts` serves as the -input for many downstream analyses. - -The input mapped reads file (`input.bam`) is in SAM/BAM format. The -reads should be sorted so those mapping to the same chromosome are -consecutive in the file. Duplicate reads should be probably be -[removed](../uniq) first, but that depends on your data. - -The methylation level for every cytosine site at single base -resolution is estimated as the ratio of methylated to total bases in -reads mapped over that site. If a site has no reads covering it, then -we leave the value at 0, but also indicate this in the output (see -below). - -Each site is marked with a *context* in the output of the `counts` -command, and you can find out more about cytosine contexts +The `counts` command takes the mapped reads and produces the methylation level +at each genomic cytosine (both strands), with the option to produce only +levels for CpG-context cytosines. While most DNA methylation in vertebrates +is in the CpG context, cytosines in other sequence contexts, such as CXG or +CHH (where H denotes adenines, thymines, or cytosines and X denotes adenines +or thymines) may also be methylated. Non-CpG methylation occurs frequently in +plant genomes and pluripotent mammalian cells such as embryonic stem +cells. And possibly whatever cells you are studying. The output of `counts` +serves as the input for many downstream analyses. + +The input mapped reads file (`input.bam`) is in SAM/BAM format. The reads +should be sorted so those mapping to the same chromosome are consecutive in +the file. Duplicate reads should be probably be removed using the +[uniq](../uniq) command first, but that depends on your data. + +The methylation level for every cytosine site at single base resolution is +estimated as the ratio of methylated to total bases in reads mapped over that +site. If a site has no reads covering it, then we leave the value at 0, but +also indicate this in the output (see below). + +Each site is marked with a *context* in the output of the `counts` command, +and you can find out more about cytosine contexts [here](../cytosine_contexts). -To compute methylation levels at each cytosine site along -the genome you can use the following command: +To compute methylation levels at each cytosine site along the genome you can +use the following command: + ```console -$ dnmtools counts -c /path/to/genome.fa -o output.meth input.sam +$ dnmtools counts -c /path/to/genome.fa -o output.counts input.sam ``` -The argument `-c` gives the name of a FASTA file containing all -chromosome sequences (as of v1.2.5, a directory of separate files is -no longer supported). Importantly, the "name" line in each chromosome -FASTA file must begin with the character `>` followed immediately by -the same name that identifies that chromosome in the SAM output (the -`.bam` files). If you use the same FASTA format file you used to map -the reads, everything should be fine. An example of the output and -explanation of each column follows: +The argument `-c` gives the name of a FASTA file containing all chromosome +sequences (as of v1.2.5, a directory of separate files is no longer +supported). Importantly, the "name" line in each chromosome FASTA file must +begin with the character `>` followed immediately by the same name that +identifies that chromosome in the SAM output (the `.bam` files). If you use +the same FASTA format file you used to map the reads, everything should be +fine. An example of the output and explanation of each column follows: + ```txt chr1 1869 + CCG 0 1 chr1 1870 + CpG 0 1 @@ -70,102 +69,96 @@ chr1 1895 - CpG 0.25 4 chr1 1898 + CHH 0 1 ``` -The output file contains one line per cytosine site. The first column -is the chromosome. The second is the position of the cytosine -(zero-based). The 3rd column indicates the strand, which can be either -'+' or '-', corresponding to either a C on the positive reference -strand or a G on the positive reference strand. The 4th column is the -sequence context of that site, followed by an x if the site has -mutated in the sample away from the reference genome. The estimate for -mutation is based on whether the reads mapping on the opposite strand -for a cytosine are biased towards 'A', which happens if a 'C' to 'T' -mutation has occurred (a G to A on the other strand) and would -otherwise be interpreted as a C-to-T conversion by bisulfite. The 5th -column is the estimated methylation level, equal to the number of Cs -in reads at position corresponding to the site, divided by the sum of -the Cs and Ts mapping to that position. The final column is number of -reads overlapping with that site. - -Note that because `counts` produces a file containing one line for -every cytosine in the genome, the file can get quite large. For -reference assembly mm10, the output is approximately 25GB. The `-n` -option produces methylation data for CpG context sites only, and for -mm10 this produces an output file that is approximately 1GB. It is -recommended that users allocate at least 8GB of memory when running -`counts`. - -To examine the methylation status of cytosines a particular sequence -context, you can use the `awk` command (or `grep`) to filter those -lines based on the fourth column. For example, to get all cytosines -within the CHG context, run the following: +The output file contains one line per cytosine site. The first column is the +chromosome. The second is the position of the cytosine (zero-based). The 3rd +column indicates the strand, which can be either '+' or '-', corresponding to +either a C on the positive reference strand or a G on the positive reference +strand. The 4th column is the sequence context of that site, followed by an x +if the site has mutated in the sample away from the reference genome. The +estimate for mutation is based on whether the reads mapping on the opposite +strand for a cytosine are biased towards 'A', which happens if a 'C' to 'T' +mutation has occurred (a G to A on the other strand) and would otherwise be +interpreted as a C-to-T conversion by bisulfite. The 5th column is the +estimated methylation level, equal to the number of Cs in reads at position +corresponding to the site, divided by the sum of the Cs and Ts mapping to that +position. The final column is number of reads overlapping with that site. + +Note that because `counts` produces a file containing one line for every +cytosine in the genome, the file can get quite large. For reference assembly +mm10, the output is approximately 25GB. The `-n` option produces methylation +data for CpG context sites only, and for mm10 this produces an output file +that is approximately 1GB. It is recommended that users allocate at least 8GB +of memory when running `counts`. + +To examine the methylation status of cytosines a particular sequence context, +you can use the `awk` command (or `grep`) to filter those lines based on the +fourth column. For example, to get all cytosines within the CHG context, run +the following: + ```console -$ awk '$4 == CHG' human_esc.meth > human_esc_chg.meth +$ awk '$4 == CHG' human_esc.counts > human_esc_chg.counts ``` -Our convention is to name `counts` output with all cytosines like -`*.meth`, with CHG like `*.chg.meth` and with CHH like `*.chh.meth`. + +Our convention is to name `counts` output with all cytosines like `*.counts`, +with CHG like `*.chg.counts` and with CHH like `*.chh.counts`. **Estimated "mutations" in the sample:** -In version 1.2.1 sites are flagged if there is an indication that a -mutation at a C in the sample might cause a T to appear in the reads, -and then confound the issues of mutations vs. the absence of -methylation. The 4th column in the output (as indicated above) is the -"context" string for that cytosine, and if a site is flagged as a -possible mutation, then an `x` will appear at the end of the context -string. This is an example of what you would see in the output: +In version 1.2.1 sites are flagged if there is an indication that a mutation +at a C in the sample might cause a T to appear in the reads, and then confound +the issues of mutations vs. the absence of methylation. The 4th column in the +output (as indicated above) is the "context" string for that cytosine, and if +a site is flagged as a possible mutation, then an `x` will appear at the end +of the context string. This is an example of what you would see in the output: ```txt chr1 1870 + CpGx 0 1 ``` -The data used to make this inference is on the opposite strand: if -there is a C on the positive strand, then we should observe a G at the -same position in reads mapping to the negative strand, regardless of -the bisulfite. If the reference has a C, but molecules in the sample -have mutated C->T, independent of methylation or bisulfite, then we -would observe an A in reads mapping to the same position on the -opposite strand. At present we use very poor statistical criteria for -this, so unless you believe you understand it well, it's best to -ignore. Hopefully we can do this in a more rigorous and interpretable -way in the future. +The data used to make this inference is on the opposite strand: if there is a +C on the positive strand, then we should observe a G at the same position in +reads mapping to the negative strand, regardless of the bisulfite. If the +reference has a C, but molecules in the sample have mutated C->T, independent +of methylation or bisulfite, then we would observe an A in reads mapping to +the same position on the opposite strand. At present we use very poor +statistical criteria for this, so unless you believe you understand it well, +it's best to ignore. Hopefully we can do this in a more rigorous and +interpretable way in the future. ## Options ```txt -o, -output ``` -Output file name. The default is to write output to the terminal, -which is not useful unless commands are piped. +Output file name. The default is to write output to the terminal, which is not +useful unless commands are piped. ```txt -c, -chrom ``` -Reference genome file, which must be in FASTA format. This is -required. +Reference genome file, which must be in FASTA format. This is required. ```txt -t, -threads ``` -The number of threads to use. This is only really helpful if the input -is BAM (not very helpful for SAM), and the output is to be zipped (see -`-z` below). These threads will help decompress the BAM input and will -help compress the gzip format output. If only one of these conditions -holds, using more threads can still help. Because `counts` spends most -of its computing time processing reads sequentially, there are +The number of threads to use. This is only helpful if the input is BAM or if +the output is to be zipped (see `-z` below). The threads speed up +decompressing BAM input or compressing gzip output. Because `counts` spends +most of its computing time processing reads sequentially, there are diminishing returns for specifying too many threads. ```txt -z, -zip ``` -The output should be zipped (in gzip format). This is not deduced by -the filename, but specifying this argument should be accompanied by -using a `.gz` filename suffix for the output. +The output should be zipped (in gzip format). This is not deduced by the +filename, but specifying this argument should be accompanied by using a `.gz` +filename suffix for the output. ```txt -n, -cpg-only ``` -Print only CpG context cytosines. This significantly reduces the -output size in most genomes. Note that using this option does not -merge data as symmetric CpGs. +Print only CpG context cytosines. This significantly reduces the output size +in most genomes. Note that using this option does not merge data as symmetric +CpGs. ```txt -v, -verbose diff --git a/docs/content/radmeth.md b/docs/content/radmeth.md index d46dc108..eb718bb9 100644 --- a/docs/content/radmeth.md +++ b/docs/content/radmeth.md @@ -3,37 +3,48 @@ ## Synopsis ```console -dnmtools radmeth -factor case -o output.txt design_matrix.txt proportion_table.txt +dnmtools radmeth -factor case -o output.txt design_matrix.txt data_matrix.txt ``` ## Description The differential-methylation detection method described in this section is based on the beta-binomial regression. We recommend using this method when -more than three replicates are available in each group. For rapid differential -methylation analysis the regression-based method should be run on a large -machine with many cores. Note that the actual processing time depends on the -coverage of each methylome, the number of sites analyzed, and the number of -methylomes in the dataset. - -### Generating proportion tables - -The first step in the differential methylation analysis is to assemble a -proportion table containing read proportions for all target -methylomes. Dnmtools includes a program [merge](../merge) to generate a -proportion table from the given collection of methylomes. Suppose that we want -to compare methylomes named `control-a.meth`, `control-b.meth`, -`control-c.meth` to the methylomes `case-a.meth`, `case-b.meth`, -`case-c.meth`. The proportion table can be created with the following -command: +more than three replicates are available in each group, or for a more +complicated experimental design, there are enough samples to fit the +parameters. The radmeth command can take advantage of multiple cores for +improved speed. The processing time depends on the number of sites analyzed +and the number of methylomes in the dataset. + +### Generating input tables + +The first step in the differential methylation analysis is to assemble a table +with the counts of total reads and methylated reads for each site (rows) and +each methylome (columns). Dnmtools includes a command [merge](../merge) to +generate this kind of table from the given collection of methylomes in +[counts](../counts) format. + +Features: although the explanation here focuses on individual CpG sites, the +identical method can be applied to genomic intervals as features (rows). The +common requirement is that in each row, each "column" has a count for total +reads and a count for reads indicating methylation. One prominent example is +to run the `roi` command on a set of genomic intervals, and use the output as +the basis for differential methylation analysis. + +Suppose that we want to compare methylomes named `control_a.sym`, +`control_b.sym`, `control_c.sym` to the methylomes `case_a.sym`, `case_b.sym`, +`case_c.sym`. Each of these filenames suggests that the methylation is for +symmetric CpG sites as would be obtained using the [sym](../sym) command. + +The data matrix can be created with the following command: ```console $ dnmtools merge -t -radmeth \ - control-a.meth control-b.meth control-c.meth \ - case-a.meth case-b.meth case-c.meth > proportion-table.txt + control_a.meth control_b.meth control_c.meth \ + case_a.meth case_b.meth case_c.meth > data_table.txt ``` -This is what `proportion-table.txt` looks like: +This is what `data_table.txt` looks like: ```txt control_a control_b control_c case_a case_b case_c @@ -43,17 +54,21 @@ chr1:160:+:CpG 12 8 10 5 17 4 15 14 1 chr1:309:+:CpG 1 1 1 0 17 12 12 8 2 1 19 8 ``` -As indicated in the header, this proportion table contains information about 6 +As indicated in the header, this table contains information about 6 methylomes: 3 controls and 3 cases. Each row of the table contains -information about a CpG site and a proportion of reads mapping over this site -in each methylome. For example, the first row describes a cytosine within a -CpG site located on chromosome 1 at position 108. This site is present in 9 -reads in the methylome control a and is methylated in 6 of them. Note that the -dnmtools `merge` command adds methylomes into the proportion table in the -order in which they are listed on the command line. - -Note: the header conforms to a traditional "data frame" from R, which does not -include a column name for the row names, which are the first column. +information about an individual CpG site, and two values for each column +heading: one value being the count of total reads and the second the count of +methylated reads. Note that the second value does not need to be integers, +which can happen when reads contribute fractionally, as in counts from +nanopore. For example, the first row describes a cytosine within a CpG site +located on chromosome 1 at position 108. This site is covered by 9 reads in +the methylome "control a" and 6 of those reads indicate methylation. Note that +the dnmtools `merge` command adds methylomes into the data matrix in the order +in which they are listed on the command line. + +Note: the header does not include a column name for the first column, and for +each methylome, currently there is only one column name for each pair of +columns in this format. ### Design matrix @@ -71,25 +86,25 @@ case_c 1 1 ``` The design matrix shows that samples in this dataset are associated with two -factors: base and case. The first column corresponds to the base factor and -will always be present in the design matrix. Think of it as stating that all +factors: base and case. The first column corresponds to the intercept and will +always be present in the design matrix. Think of it as stating that all samples have the same baseline mean methylation level. To distinguish cases from controls we add another factor case (second column). The `1`s in this -column correspond to the samples which belong to the cases group. You can use -this design matrix as a template to create design matrices for two-group +second column indicate the methylomes that belong to the cases group. You can +use this design matrix as a template to create design matrices for two-group comparisons involving arbitrary many methylomes. -After creating the proportion table and the design matrix, we are now ready to -start the differential methylation analysis. It consists of (1) regression, +After creating the data matrix and the design matrix, we are now ready to +start the differential methylation analysis. It consists of (1) regression, (2) combining significance, and (3) multiple testing adjustment steps. ### Regression -Suppose that the `proportion-table.txt` and `design-matrix.txt` are as +Suppose that the `data_table.txt` and `design_matrix.txt` are as described above. The regression step is run with the command ```console -$ dnmtools radmeth -factor case -o output.txt design-matrix.txt proportion-table.txt +$ dnmtools radmeth -factor case -o output.txt design_matrix.txt data_table.txt ``` The `-factor` parameter specifies the factor with respect to which we want to @@ -101,18 +116,36 @@ respectively. The p-value (5-th column) is given by the log-likelihood ratio test. Depending on the distribution of methylated and unmethylated CpGs, some -p-values may not be calculated. If there are no reads that cover a given CpG, -the value `NA_LOW_COV` will be printed instead of a p-value. If all CpGs are -fully methylated and fully methylated in both case and control, we say the CpG -has an "extreme distribution". In these cases, we cannot reject the null and -output `NA_EXTREME_CNT`. Finally, if there are numerical errors in the -log-likelihood ratio test, we simply print `NA`. You can opt out of these -additional NA informations by not using the `-n` flag. - -We do not recommend relying too much on the individual p-values of CpGs from -radmeth, as they often do not contain sufficient coverage for significant -biological interpretations. Instead, we recommend merging p-values of -neighboring CpGs using [radadjust](../radadjust). +p-values may not be calculated, and a value of NA will be output. If the +argument `-na-info` is used, the following detail will be given. If there are +no reads that cover a given CpG, the value `NA_LOW_COV` will be printed +instead of a p-value. If all CpGs are fully methylated and fully methylated in +both case and control, we say the CpG has a "extreme" counts. In these cases, +we cannot reject the null and output `NA_EXTREME_CNT`. Finally, if there are +numerical errors in the log-likelihood ratio test, we simply print `NA`. + +## Output format + +The output for `radmeth` is one line for each row of the input data matrix. +If the input data matrix has rownames with "tokens" separated by colons, for +example "chr3:67363494:+:CpG", then this will appear in the output but with +the colons replaced by tabs. If the rownames in the input file do not have +colons, they will be output as is, in the first column of the output file. + +- The second column, or the first column after the (possibly transformed) + rowname, is the p-value for that column, or an NA value (see above). +- After the p-value, the mean methylation per "group" is given in each + subsequent column. Here, a "group" is a combination of factor levels, + assuming binary factors. The order of groups corresponds to the order of + their binary representation in the design matrix. The order for the groups + in the example design matrix below would be (MB=100, MA=101, FB=110, + FA=111). +- The final column in the output is the "overdispersion factor" for the + beta-binomial. Assume `phi` is the dispersion. The binomial variance is + `np(1-p)` and the beta-binomial variance is `np(1-p)((n+phi)/(phi+1))`. So + `(n - 1)/(phi + 1)` is the factor by which the variance for the + beta-binomial exceeds that of the binomial. An overdispersion factor of 0 + means the beta-binomial is behaving just like a binomial. ## Tutorial @@ -281,17 +314,17 @@ differential methylation between levels for our factor of interest (named "factor" in the design matrix). Here is a look at part of the output: ```console -$ head sim/samples.radmeth -chr1 163 + CG 0.0877822 35 22 45 36 -chr1 206 + CG 0.442745 37 30 43 38 -chr1 232 + CG 0.756925 28 20 35 26 -chr1 278 + CG 0.794626 43 35 37 31 -chr1 296 + CG 0.635986 33 29 38 32 -chr1 310 + CG 0.641986 37 28 31 22 -chr1 322 + CG 0.0973903 38 33 47 35 -chr1 324 + CG 0.260703 40 29 31 26 -chr1 350 + CG 0.888242 36 30 35 29 -chr1 356 + CG 0.0535286 37 35 42 33 +$ head samples.radmeth | column -t +chr1 163 + CG 0.0863989 0.818157 0.653660 0.775035 0.591029 0.000052 +chr1 206 + CG 0.443989 0.841557 0.764998 0.908700 0.859154 0.000091 +chr1 232 + CG 0.756234 0.776885 0.744420 0.725115 0.688140 0.000242 +chr1 278 + CG 0.77863 0.831485 0.807046 0.842138 0.818908 0.000808 +chr1 296 + CG 0.648263 0.848138 0.884828 0.836674 0.875727 0.000341 +chr1 310 + CG 0.640836 0.729179 0.780351 0.689915 0.745921 0.139695 +chr1 322 + CG 0.0972604 0.679385 0.849711 0.789027 0.908915 0.000235 +chr1 324 + CG 0.25803 0.861294 0.753022 0.822346 0.694456 0.008238 +chr1 350 + CG 0.890771 0.844111 0.855183 0.807851 0.820951 0.000267 +chr1 356 + CG 0.0527503 0.750013 0.933344 0.818157 0.954543 0.000315 ``` Focusing on the 5th column, we see values between 0 and 1, but none of them @@ -306,16 +339,19 @@ for i in M F A B; do done ``` +(The argument `-relaxed` to the [selectsites](../selectsites) command is to +allow for flexibility in the column format after the first 6.) + The result allows us to verify that `radmeth` with the design matrix given above identifies sites as significant when they differ between levels A and B of the factor of interest: ```console $ for i in sim/features_[ABFM].txt; do echo $i `awk 'BEGIN{k=0;c=0}{k+=$5;c+=1}END{print k,c,k/c}' $i`; done | column -t -sim/features_A.txt 0.256072 430 0.000595517 -sim/features_B.txt 0.114786 372 0.000308564 -sim/features_F.txt 197.332 400 0.493329 -sim/features_M.txt 297.728 588 0.50634 +sim/features_A.txt 0.255514 430 0.000594218 +sim/features_B.txt 0.11431 372 0.000307286 +sim/features_F.txt 196.602 400 0.491505 +sim/features_M.txt 296.222 588 0.50378 ``` The p-values in features associated with either M or F are averaging around @@ -354,10 +390,10 @@ analogous but with this modified design matrix: $ dnmtools radmeth -o sim/samples_noint.radmeth -f factor design_noint.txt sim/table.txt $ for i in M F A B; do dnmtools selectsites -relaxed -o sim/features_noint_${i}.txt sim/features_${i}.bed sim/samples_noint.radmeth; done $ for i in sim/*_noint_*.txt; do echo $i `awk 'BEGIN{k=0;c=0}{k+=$5;c+=1}END{print k,c,k/c}' $i`; done | column -t -sim/features_noint_A.txt 8.27931 430 0.0192542 -sim/features_noint_B.txt 15.2397 372 0.040967 -sim/features_noint_F.txt 77.2656 400 0.193164 -sim/features_noint_M.txt 93.7312 588 0.159407 +sim/features_noint_A.txt 8.27797 430 0.0192511 +sim/features_noint_B.txt 15.2391 372 0.0409654 +sim/features_noint_F.txt 77.2454 400 0.193114 +sim/features_noint_M.txt 93.7184 588 0.159385 ``` Although the sites within features corresponding to A and B are much lower on diff --git a/docs/content/selectsites.md b/docs/content/selectsites.md index 2e812109..6041b2f3 100644 --- a/docs/content/selectsites.md +++ b/docs/content/selectsites.md @@ -1,8 +1,8 @@ -# selectsites - Subsetting cytosines in methcounts files +# selectsites - get subsets of cytosines from counts files ## Synopsis ```shell - $ dnmtools selectsites [OPTIONS] +$ dnmtools selectsites [OPTIONS] ``` ## Description @@ -13,19 +13,17 @@ calculating average methylation levels in (1) annotated regions, such as promoter regions or repeats or (2) regions defined by the data itself, such as HMRs or PMDs. -A possible solution to subset these regions is to convert the -methcounts file to BED format, intersect it with a BED file of the -regions of interest (using +A possible solution to subset these regions is to convert the counts file to +BED format, intersect it with a BED file of the regions of interest (using [bedtools](https://bedtools.readthedocs.io)), then convert it back to -methcounts. The program selectsites simplifies these operations. It -takes a [counts](../counts) file and a set of intervals -represented as a BED file and produces a subset of the entries in the -methcounts file included in the BED regions. We can select entries in -input.meth contained in any inverval in `regions.bed` using the -following command. +counts. The program selectsites simplifies these operations. It takes a +[counts](../counts) format file and a set of intervals in a BED file and +produces a subset of the entries in the counts file included in the BED +regions. We can select entries in `input.counts` contained in any inverval in +`regions.bed` using the following command. ```shell -$ dnmtools selectsites -o output.meth regions.bed input.meth +$ dnmtools selectsites -o output.counts regions.bed input.counts ``` ## Options @@ -38,24 +36,29 @@ Name of output file (default: STDOUT) ```txt -p, -preload ``` -preload sites (use for large target intervals) +Preload sites (use for large target intervals). ```txt -v, -verbose ``` -print more run info to STDERR while the program is running. +Print more run info to STDERR while the program is running. ```txt -d, -disk ``` -process sites on disk (fast if target intervals are few) +Process sites on disk (fast if target intervals are few). ```txt -S, -summary ``` -write summary to this file +Write summary to this file. ```txt -z, -zip ``` -output file will be in gzip compressed format +The output file will be in gzip compressed format. + +```txt + -relaxed +``` +Allow additional columns in the input file. diff --git a/src/analysis/methcounts.cpp b/src/analysis/methcounts.cpp index 469d44c1..9caa219a 100644 --- a/src/analysis/methcounts.cpp +++ b/src/analysis/methcounts.cpp @@ -1,22 +1,30 @@ -/* counts: a program for counting the methylated and unmethylated - * reads mapping over each CpG or C +/* counts: a program for counting the methylated and unmethylated reads + * mapping over each CpG or C * - * Copyright (C) 2011-2023 University of Southern California and - * Andrew D. Smith + * Copyright (C) 2011-2025 Andrew D. Smith * * Authors: Andrew D. Smith, Song Qiang, Guilherme Sena, and Masaru Nakajima * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. */ +#include "bam_record_utils.hpp" +#include "bsutils.hpp" +#include "counts_header.hpp" + +#include "OptionParser.hpp" + +/* HTSlib */ +#include + #include // for [u]int[0-9]+_t #include #include @@ -26,66 +34,20 @@ #include #include -#include "OptionParser.hpp" -// #include "GenomicRegion.hpp" -/* ADS: This code writes MSite objects to files, but does not use - MSite to do it. Possiby dangerous, but currently much faster. If - MSite has a way to serialize into a char[] directly, then we should - use it. */ -// #include "MSite.hpp" -#include "bam_record_utils.hpp" -#include "bsutils.hpp" -#include "counts_header.hpp" -#include "dnmt_error.hpp" - -/* HTSlib */ -#include - -using std::cerr; -using std::endl; -using std::string; -using std::unordered_map; -using std::unordered_set; -using std::vector; - -using bamxx::bam_header; -using bamxx::bam_rec; -using bamxx::bgzf_file; - -struct quick_buf : public std::ostringstream, - public std::basic_stringbuf { - // ADS: By user ecatmur on SO; very fast. Seems to work... - quick_buf() { - // ...but this seems to depend on data layout - static_cast &>(*this).rdbuf(this); - } - void - clear() { - // reset buffer pointers (member functions) - setp(pbase(), pbase()); - } - char const * - c_str() { - /* between c_str and insertion make sure to clear() */ - *pptr() = '\0'; - return pbase(); - } -}; - // ADS: we should never have to worry about coverage over > 32767 in // any downstream analysis, so using "int16_t" here would allow to // detect wrap around and report it as some kind of weird thing, maybe // zeroing it and flagging the output. As it is now, if we really need // >65535-fold coverage, we can make the change here. -typedef uint16_t count_type; +typedef std::uint16_t count_type; -static inline bool -eats_ref(const uint32_t c) { +[[nodiscard]] static inline bool +eats_ref(const std::uint32_t c) { return bam_cigar_type(bam_cigar_op(c)) & 2; } -static inline bool -eats_query(const uint32_t c) { +[[nodiscard]] static inline bool +eats_query(const std::uint32_t c) { return bam_cigar_type(bam_cigar_op(c)) & 1; } @@ -95,33 +57,31 @@ eats_query(const uint32_t c) { then one needs this. Also, Qiang should be consulted on this because he spent much time thinking about it in the context of plants. */ -static bool +[[nodiscard]] static bool is_chh(const std::string &s, size_t i) { return (i < (s.length() - 2)) && is_cytosine(s[i]) && !is_guanine(s[i + 1]) && !is_guanine(s[i + 2]); } -static bool +[[nodiscard]] static bool is_ddg(const std::string &s, size_t i) { return (i < (s.length() - 2)) && !is_cytosine(s[i]) && !is_cytosine(s[i + 1]) && is_guanine(s[i + 2]); } -static bool +[[nodiscard]] static bool is_c_at_g(const std::string &s, size_t i) { return (i < (s.length() - 2)) && is_cytosine(s[i]) && !is_cytosine(s[i + 1]) && !is_guanine(s[i + 1]) && is_guanine(s[i + 2]); } -/* Right now the CountSet objects below are much larger than they need - to be, for the things we are computing. However, it's not clear - that the minimum information would really put the memory - requirement of the program into a more reasonable range, so keeping - all the information seems reasonable. */ +/* Right now the CountSet objects below are much larger than they need to be, + for the things we are computing. However, it's not clear that the minimum + information would really put the memory requirement of the program into a + more reasonable range, so keeping all the information seems reasonable. */ struct CountSet { - - string + std::string tostring() const { std::ostringstream oss; oss << pA << '\t' << pC << '\t' << pG << '\t' << pT << '\t' << nA << '\t' @@ -184,14 +144,13 @@ struct CountSet { // count_type N; /* this wasn't used and breaks alignment */ }; -/* The "tag" returned by this function should be exclusive, so that - * the order of checking conditions doesn't matter. There is also a - * bit of a hack in that the unsigned "pos" could wrap, but this still - * works as long as the chromosome size is not the maximum size of a - * size_t. +/* The "tag" returned by this function should be exclusive, so that the order + * of checking conditions doesn't matter. There is also a bit of a hack in + * that the unsigned "pos" could wrap, but this still works as long as the + * chromosome size is not the maximum size of a size_t. */ -static uint32_t -get_tag_from_genome(const string &s, const size_t pos) { +[[nodiscard]] static std::uint32_t +get_tag_from_genome(const std::string &s, const size_t pos) { if (is_cytosine(s[pos])) { if (is_cpg(s, pos)) return 0; @@ -215,16 +174,16 @@ get_tag_from_genome(const string &s, const size_t pos) { return 4; // shouldn't be used for anything } -/* This "has_mutated" function looks on the opposite strand to see - * if the apparent conversion from C->T was actually already in the - * DNA because of a mutation or SNP. +/* This "has_mutated" function looks on the opposite strand to see if the + * apparent conversion from C->T was actually already in the DNA because of a + * mutation or SNP. */ -static bool +[[nodiscard]] static bool has_mutated(const char base, const CountSet &cs) { - static const double MUTATION_DEFINING_FRACTION = 0.5; + static constexpr auto mutation_defining_frac = 0.5; return is_cytosine(base) - ? (cs.nG < MUTATION_DEFINING_FRACTION * (cs.neg_total())) - : (cs.pG < MUTATION_DEFINING_FRACTION * (cs.pos_total())); + ? (cs.nG < mutation_defining_frac * (cs.neg_total())) + : (cs.pG < mutation_defining_frac * (cs.pos_total())); } static const char *tag_values[] = { @@ -239,26 +198,27 @@ static const char *tag_values[] = { "CCGx", // 8 "Nx" // 9 }; -static const uint32_t MUT_OFFSET = 5; +static const std::uint32_t MUT_OFFSET = 5; -static inline uint32_t -tag_with_mut(const uint32_t tag, const bool mut) { +[[nodiscard]] static inline std::uint32_t +tag_with_mut(const std::uint32_t tag, const bool mut) { return tag + (mut ? MUT_OFFSET : 0); } template static void -write_output(const bam_header &hdr, bgzf_file &out, const int32_t tid, - const string &chrom, const vector &counts, - bool CPG_ONLY) { - - quick_buf buf; // keep underlying buffer space? - - for (size_t i = 0; i < counts.size(); ++i) { +write_output(const bamxx::bam_header &hdr, bamxx::bgzf_file &out, + const int32_t tid, const std::string &chrom, + const std::vector &counts, bool CPG_ONLY) { + constexpr auto buf_size = 1024; // max width of a line for counts output + constexpr const char *fmt = "%s\t%ld\t%c\t%s\t%.6g\t%d\n"; + char buf[buf_size]; + + for (size_t i = 0; i < size(counts); ++i) { const char base = chrom[i]; if (is_cytosine(base) || is_guanine(base)) { - const uint32_t the_tag = get_tag_from_genome(chrom, i); + const std::uint32_t the_tag = get_tag_from_genome(chrom, i); if (CPG_ONLY && the_tag != 0) continue; @@ -267,27 +227,29 @@ write_output(const bam_header &hdr, bgzf_file &out, const int32_t tid, : counts[i].unconverted_guanine(); const double converted = is_c ? counts[i].converted_cytosine() : counts[i].converted_guanine(); - const uint32_t n_reads = unconverted + converted; + const std::uint32_t n_reads = unconverted + converted; if (require_covered && n_reads == 0) continue; const bool mut = has_mutated(base, counts[i]); - buf.clear(); - // ADS: here is where we make an MSite, but not using MSite - buf << sam_hdr_tid2name(hdr, tid) << '\t' << i << '\t' - << (is_c ? '+' : '-') << '\t' - << tag_values[tag_with_mut(the_tag, mut)] << '\t' - << (n_reads > 0 ? unconverted / n_reads : 0.0) << '\t' << n_reads - << '\n'; - if (!out.write(buf.c_str(), buf.tellp())) - throw dnmt_error("error writing output"); + // clang-format off + const int n = std::snprintf(buf, buf_size, fmt, + sam_hdr_tid2name_ptr(hdr, tid), + i, + (is_c ? '+' : '-'), + tag_values[tag_with_mut(the_tag, mut)], + (n_reads > 0 ? unconverted / n_reads : 0.0), + n_reads); + // clang-format on + if (n < 0 || !out.write(buf, n)) + throw std::runtime_error("error formatting output"); } } } static void -count_states_pos(const bam_rec &aln, vector &counts) { +count_states_pos(const bamxx::bam_rec &aln, std::vector &counts) { /* Move through cigar, reference and read positions without inflating cigar or read sequence */ const auto seq = bam_get_seq(aln); @@ -297,7 +259,7 @@ count_states_pos(const bam_rec &aln, vector &counts) { auto qpos = 0; // to match type with b->core.l_qseq for (auto c_itr = beg_cig; c_itr != end_cig; ++c_itr) { const char op = bam_cigar_op(*c_itr); - const uint32_t n = bam_cigar_oplen(*c_itr); + const std::uint32_t n = bam_cigar_oplen(*c_itr); if (eats_ref(op) && eats_query(op)) { const decltype(qpos) end_qpos = qpos + n; for (; qpos < end_qpos; ++qpos) { @@ -314,15 +276,15 @@ count_states_pos(const bam_rec &aln, vector &counts) { rpos += n; } } - // ADS: somehow previous code included a correction for rpos going - // past the end of the chromosome; this should result at least in a - // soft-clip by any mapper. I'm not checking it here as even if it - // happens I don't want to terminate. + // ADS: somehow previous code included a correction for rpos going past the + // end of the chromosome; this should result at least in a soft-clip by any + // mapper. I'm not checking it here as even if it happens I don't want to + // terminate. assert(qpos == get_l_qseq(aln)); } static void -count_states_neg(const bam_rec &aln, vector &counts) { +count_states_neg(const bamxx::bam_rec &aln, std::vector &counts) { /* Move through cigar, reference and (*backward*) through read positions without inflating cigar or read sequence */ const auto seq = bam_get_seq(aln); @@ -332,7 +294,7 @@ count_states_neg(const bam_rec &aln, vector &counts) { size_t qpos = get_l_qseq(aln); // to match type with b->core.l_qseq for (auto c_itr = beg_cig; c_itr != end_cig; ++c_itr) { const char op = bam_cigar_op(*c_itr); - const uint32_t n = bam_cigar_oplen(*c_itr); + const std::uint32_t n = bam_cigar_oplen(*c_itr); if (eats_ref(op) && eats_query(op)) { const size_t end_qpos = qpos - n; // to match type with qpos for (; qpos > end_qpos; --qpos) // beware ++ in macro below!!! @@ -350,17 +312,17 @@ count_states_neg(const bam_rec &aln, vector &counts) { assert(qpos == 0); } -static unordered_map -get_tid_to_idx(const bam_header &hdr, - const unordered_map name_to_idx) { - unordered_map tid_to_idx; +[[nodiscard]] static std::unordered_map +get_tid_to_idx(const bamxx::bam_header &hdr, + const std::unordered_map name_to_idx) { + std::unordered_map tid_to_idx; for (int32_t i = 0; i < hdr.h->n_targets; ++i) { // "curr_name" gives a "tid_to_name" mapping allowing to jump // through "name_to_idx" and get "tid_to_idx" - const string curr_name(hdr.h->target_name[i]); + const std::string curr_name(hdr.h->target_name[i]); const auto name_itr(name_to_idx.find(curr_name)); - if (name_itr == end(name_to_idx)) - throw dnmt_error("failed to find chrom: " + curr_name); + if (name_itr == std::cend(name_to_idx)) + throw std::runtime_error("failed to find chrom: " + curr_name); tid_to_idx[i] = name_itr->second; } return tid_to_idx; @@ -368,17 +330,18 @@ get_tid_to_idx(const bam_header &hdr, template static void -output_skipped_chromosome(const bool CPG_ONLY, const int32_t tid, - const unordered_map &tid_to_idx, - const bam_header &hdr, - const vector::const_iterator chroms_beg, - const vector &chrom_sizes, vector &counts, - bgzf_file &out) { +output_skipped_chromosome( + const bool CPG_ONLY, const int32_t tid, + const std::unordered_map &tid_to_idx, + const bamxx::bam_header &hdr, + const std::vector::const_iterator chroms_beg, + const std::vector &chrom_sizes, std::vector &counts, + bamxx::bgzf_file &out) { // get the index of the next chrom sequence const auto chrom_idx = tid_to_idx.find(tid); - if (chrom_idx == cend(tid_to_idx)) - throw dnmt_error("chrom not found: " + sam_hdr_tid2name(hdr, tid)); + if (chrom_idx == std::cend(tid_to_idx)) + throw std::runtime_error("chrom not found: " + sam_hdr_tid2name(hdr, tid)); const auto chrom_itr = chroms_beg + chrom_idx->second; @@ -389,19 +352,20 @@ output_skipped_chromosome(const bool CPG_ONLY, const int32_t tid, write_output(hdr, out, tid, *chrom_itr, counts, CPG_ONLY); } -static bool -consistent_targets(const bam_header &hdr, - const unordered_map &tid_to_idx, - const vector &names, const vector &sizes) { +[[nodiscard]] static bool +consistent_targets(const bamxx::bam_header &hdr, + const std::unordered_map &tid_to_idx, + const std::vector &names, + const std::vector &sizes) { const size_t n_targets = hdr.h->n_targets; if (n_targets != names.size()) return false; for (size_t tid = 0; tid < n_targets; ++tid) { - const string tid_name_sam = sam_hdr_tid2name(hdr, tid); + const std::string tid_name_sam = sam_hdr_tid2name(hdr, tid); const size_t tid_size_sam = sam_hdr_tid2len(hdr, tid); const auto idx_itr = tid_to_idx.find(tid); - if (idx_itr == cend(tid_to_idx)) + if (idx_itr == std::cend(tid_to_idx)) return false; const auto idx = idx_itr->second; if (tid_name_sam != names[idx] || tid_size_sam != sizes[idx]) @@ -414,21 +378,21 @@ template static void process_reads(const bool VERBOSE, const bool show_progress, const bool compress_output, const bool include_header, - const size_t n_threads, const string &infile, - const string &outfile, const string &chroms_file, + const size_t n_threads, const std::string &infile, + const std::string &outfile, const std::string &chroms_file, const bool CPG_ONLY) { // first get the chromosome names and sequences from the FASTA file - vector chroms, names; + std::vector chroms, names; read_fasta_file_short_names(chroms_file, names, chroms); for (auto &i : chroms) - transform(cbegin(i), cend(i), begin(i), - [](const char c) { return std::toupper(c); }); + std::transform(std::cbegin(i), std::cend(i), std::begin(i), + [](const char c) { return std::toupper(c); }); if (VERBOSE) - cerr << "[n chroms in reference: " << chroms.size() << "]" << endl; - const auto chroms_beg = cbegin(chroms); + std::cerr << "[n chroms in reference: " << chroms.size() << "]\n"; + const auto chroms_beg = std::cbegin(chroms); - unordered_map name_to_idx; - vector chrom_sizes(chroms.size(), 0); + std::unordered_map name_to_idx; + std::vector chrom_sizes(chroms.size(), 0); for (size_t i = 0; i < chroms.size(); ++i) { name_to_idx[names[i]] = i; chrom_sizes[i] = chroms[i].size(); @@ -439,22 +403,23 @@ process_reads(const bool VERBOSE, const bool show_progress, // open the hts SAM/BAM input file and get the header bamxx::bam_in hts(infile); if (!hts) - throw dnmt_error("failed to open input file"); + throw std::runtime_error("failed to open input file"); // load the input file's header - bam_header hdr(hts); + bamxx::bam_header hdr(hts); if (!hdr) - throw dnmt_error("failed to read header"); + throw std::runtime_error("failed to read header"); - unordered_map tid_to_idx = get_tid_to_idx(hdr, name_to_idx); + std::unordered_map tid_to_idx = + get_tid_to_idx(hdr, name_to_idx); if (!consistent_targets(hdr, tid_to_idx, names, chrom_sizes)) - throw dnmt_error("inconsistent reference genome information"); + throw std::runtime_error("inconsistent reference genome information"); // open the output file - const string output_mode = compress_output ? "w" : "wu"; - bgzf_file out(outfile, output_mode); + const std::string output_mode = compress_output ? "w" : "wu"; + bamxx::bgzf_file out(outfile, output_mode); if (!out) - throw dnmt_error("error opening output file: " + outfile); + throw std::runtime_error("error opening output file: " + outfile); // set the threads for the input file decompression if (n_threads > 1) { @@ -467,13 +432,13 @@ process_reads(const bool VERBOSE, const bool show_progress, // now iterate over the reads, switching chromosomes and writing // output as needed - bam_rec aln; + bamxx::bam_rec aln; int32_t prev_tid = -1; // this is where all the counts are accumulated - vector counts; + std::vector counts; - vector::const_iterator chrom_itr; + std::vector::const_iterator chrom_itr; while (hts.read(hdr, aln)) { const int32_t tid = get_tid(aln); @@ -492,7 +457,7 @@ process_reads(const bool VERBOSE, const bool show_progress, "previous tid: " + std::to_string(prev_tid) + " current tid: " + std::to_string(tid); - throw dnmt_error(message); + throw std::runtime_error(message); } if (!require_covered) @@ -502,11 +467,11 @@ process_reads(const bool VERBOSE, const bool show_progress, // get the next chrom to process auto chrom_idx(tid_to_idx.find(tid)); - if (chrom_idx == end(tid_to_idx)) - throw dnmt_error("chromosome not found: " + - string(sam_hdr_tid2name(hdr, tid))); + if (chrom_idx == std::cend(tid_to_idx)) + throw std::runtime_error("chromosome not found: " + + std::string(sam_hdr_tid2name(hdr, tid))); if (show_progress) - cerr << "processing " << sam_hdr_tid2name(hdr, tid) << endl; + std::cerr << "processing " << sam_hdr_tid2name(hdr, tid) << '\n'; prev_tid = tid; chrom_itr = chroms_beg + chrom_idx->second; @@ -525,8 +490,8 @@ process_reads(const bool VERBOSE, const bool show_progress, write_output(hdr, out, prev_tid, *chrom_itr, counts, CPG_ONLY); - // ADS: if some chroms might not be covered by reads, we have to - // iterate over what remains + // ADS: if some chroms might not be covered by reads, we have to iterate + // over what remains if (!require_covered) for (auto i = prev_tid + 1; i < hdr.h->n_targets; ++i) output_skipped_chromosome(CPG_ONLY, i, tid_to_idx, hdr, chroms_beg, @@ -545,8 +510,8 @@ main_counts(int argc, char *argv[]) { bool compress_output = false; bool include_header = false; - string chroms_file; - string outfile; + std::string chroms_file; + std::string outfile; int n_threads = 1; /****************** COMMAND LINE OPTIONS ********************/ @@ -569,40 +534,40 @@ main_counts(int argc, char *argv[]) { opt_parse.add_opt("zip", 'z', "output gzip format", false, compress_output); opt_parse.add_opt("progress", '\0', "show progress", false, show_progress); opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE); - vector leftover_args; + std::vector leftover_args; opt_parse.parse(argc, argv, leftover_args); if (opt_parse.about_requested() || opt_parse.help_requested() || leftover_args.empty()) { - cerr << opt_parse.help_message() << endl - << opt_parse.about_message() << endl; + std::cerr << opt_parse.help_message() << '\n' + << opt_parse.about_message() << '\n'; return EXIT_SUCCESS; } if (opt_parse.option_missing()) { - cerr << opt_parse.option_missing_message() << endl; + std::cerr << opt_parse.option_missing_message() << '\n'; return EXIT_SUCCESS; } - const string mapped_reads_file = leftover_args.front(); + const std::string mapped_reads_file = leftover_args.front(); /****************** END COMMAND LINE OPTIONS *****************/ if (n_threads < 0) - throw dnmt_error("thread count cannot be negative"); + throw std::runtime_error("thread count cannot be negative"); std::ostringstream cmd; - copy(argv, argv + argc, std::ostream_iterator(cmd, " ")); + std::copy(argv, argv + argc, std::ostream_iterator(cmd, " ")); // file types from HTSlib use "-" for the filename to go to stdout if (outfile.empty()) outfile = "-"; if (VERBOSE) - cerr << "[input BAM/SAM file: " << mapped_reads_file << "]" << endl - << "[output file: " << outfile << "]" << endl - << "[output format: " << (compress_output ? "bgzf" : "text") << "]" - << endl - << "[genome file: " << chroms_file << "]" << endl - << "[threads requested: " << n_threads << "]" << endl - << "[CpG only mode: " << (CPG_ONLY ? "yes" : "no") << "]" << endl - << "[command line: \"" << cmd.str() << "\"]" << endl; + std::cerr << "[input BAM/SAM file: " << mapped_reads_file << "]\n" + << "[output file: " << outfile << "]\n" + << "[output format: " << (compress_output ? "bgzf" : "text") + << "]\n" + << "[genome file: " << chroms_file << "]\n" + << "[threads requested: " << n_threads << "]\n" + << "[CpG only mode: " << (CPG_ONLY ? "yes" : "no") << "]\n" + << "[command line: \"" << cmd.str() << "\"]\n"; if (require_covered) process_reads(VERBOSE, show_progress, compress_output, @@ -614,7 +579,7 @@ main_counts(int argc, char *argv[]) { CPG_ONLY); } catch (const std::exception &e) { - cerr << e.what() << endl; + std::cerr << e.what() << '\n'; return EXIT_FAILURE; } return EXIT_SUCCESS; diff --git a/src/common/MSite.cpp b/src/common/MSite.cpp index 14a67131..d4abd2fc 100644 --- a/src/common/MSite.cpp +++ b/src/common/MSite.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include "counts_header.hpp" #include "smithlab_utils.hpp" @@ -50,7 +51,7 @@ MSite::initialize(const char *c, const char *c_end) { failed = true; { - const uint32_t d = std::distance(field_s, field_e); + const std::uint32_t d = std::distance(field_s, field_e); chrom = string{field_s, d}; } @@ -75,7 +76,7 @@ MSite::initialize(const char *c, const char *c_end) { failed = failed || (field_e == c_end); { - const uint32_t d = std::distance(field_s, field_e); + const std::uint32_t d = std::distance(field_s, field_e); context = string{field_s, d}; } @@ -217,12 +218,10 @@ find_offset_for_msite(const std::string &chr, const size_t idx, } } -#include -using std::unordered_map; void -find_offset_for_msite(const unordered_map &chrom_order, - const std::string &chr, const size_t idx, - std::ifstream &site_in) { +find_offset_for_msite( + const std::unordered_map &chrom_order, + const std::string &chr, const size_t idx, std::ifstream &site_in) { site_in.seekg(0, ios_base::beg); const size_t begin_pos = site_in.tellg(); diff --git a/src/common/bam_record_utils.hpp b/src/common/bam_record_utils.hpp index 7fbd4a42..46221960 100644 --- a/src/common/bam_record_utils.hpp +++ b/src/common/bam_record_utils.hpp @@ -274,6 +274,11 @@ sam_hdr_tid2name(const bamxx::bam_header &hdr, const int32_t tid) { return std::string(sam_hdr_tid2name(hdr.h, tid)); } +inline const char * +sam_hdr_tid2name_ptr(const bamxx::bam_header &hdr, const int32_t tid) { + return sam_hdr_tid2name(hdr.h, tid); +} + inline uint32_t sam_hdr_tid2len(const bamxx::bam_header &hdr, const int32_t tid) { return sam_hdr_tid2len(hdr.h, tid); diff --git a/src/radmeth/radmeth.cpp b/src/radmeth/radmeth.cpp index 789304c5..fcf9bae1 100644 --- a/src/radmeth/radmeth.cpp +++ b/src/radmeth/radmeth.cpp @@ -3,15 +3,15 @@ * Author: Andrew D. Smith * Contributors: Egor Dolzhenko and Guilherme Sena * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. + * This program is free software: you can redistribute it and/or modify it + * under the terms of the GNU General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) + * any later version. * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. + * This program is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for + * more details. */ #include "radmeth_model.hpp" @@ -242,6 +242,11 @@ read_design(const std::string &design_filename) { return design; } +[[nodiscard]] static inline double +overdispersion_factor(const std::uint32_t n_samples, const double dispersion) { + return (n_samples - 1) / (dispersion + 1); +} + enum class row_status : std::uint8_t { ok, na, @@ -327,7 +332,7 @@ that the design matrix and the proportion table are correctly formatted. if (t_alt_model.props_size() != n_samples) throw std::runtime_error("found row with wrong number of columns"); - const auto [p_val, status] = [&]() -> std::tuple { + const auto p_val_status = [&]() -> std::tuple { // Skip the test if (1) no coverage in all cases or in all // controls, or (2) the site is completely methylated or // completely unmethylated across all samples. @@ -349,12 +354,14 @@ that the design matrix and the proportion table are correctly formatted. return (p_value != p_value) ? std::tuple{1.0, row_status::na} : std::tuple{p_value, row_status::ok}; }(); + // ADS: avoid capture structured binding in C++17 + const auto p_val = std::get<0>(p_val_status); + const auto status = std::get<1>(p_val_status); n_bytes[b] = [&] { // clang-format off const int n_prefix_bytes = - std::sprintf(bufs[b].data(), "%s\t", - t_alt_model.props.rowname.data()); + std::sprintf(bufs[b].data(), "%s\t", t_alt_model.rowname().data()); // clang-format on if (n_prefix_bytes < 0) return n_prefix_bytes; @@ -388,7 +395,8 @@ that the design matrix and the proportion table are correctly formatted. return n; n_bytes += n; } - const int n = std::sprintf(cursor, "\t%f\n", phi_estim_alt); + const auto od = overdispersion_factor(n_samples, phi_estim_alt); + const int n = std::sprintf(cursor, "\t%f\n", od); if (n < 0) return n; n_bytes += n; @@ -534,7 +542,7 @@ main_radmeth(int argc, char *argv[]) { << "(5) p-value\n" << "(6) estimated methylation 0\n" << "(7) estimated methylation 1\n" - << "(8) estimated dispersion\n" + << "(8) overdispersion factor (variance above binomial)\n" << "estimated methylation is for test factor value (0 or 1)\n" << '\n'; // clang-format on diff --git a/src/radmeth/radmeth_model.cpp b/src/radmeth/radmeth_model.cpp index 0df305cc..1aec01fe 100644 --- a/src/radmeth/radmeth_model.cpp +++ b/src/radmeth/radmeth_model.cpp @@ -56,14 +56,20 @@ SiteProp::parse(const std::string &line) { field_e = std::find_if(field_s + 1, line_end, is_sep); { const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_reads); - failed = failed || (ptr != field_e); + failed = failed || (ec != std::errc{}); } + // get the methylated count field_s = std::find_if(field_e + 1, line_end, not_sep); field_e = std::find_if(field_s + 1, line_end, is_sep); { +#ifdef __APPLE__ + const int ret = std::sscanf(field_s, "%lf", &mc1.n_meth); + failed = failed || (ret < 1); +#else const auto [ptr, ec] = std::from_chars(field_s, field_e, mc1.n_meth); - failed = failed || (ptr != field_e); + failed = failed || (ec != std::errc{}); +#endif } mc.push_back(mc1); diff --git a/src/radmeth/radmeth_model.hpp b/src/radmeth/radmeth_model.hpp index fb4b5ac9..013e51f4 100644 --- a/src/radmeth/radmeth_model.hpp +++ b/src/radmeth/radmeth_model.hpp @@ -25,12 +25,6 @@ #include #include -struct cumul_counts { - std::vector m_counts; - std::vector r_counts; - std::vector d_counts; -}; - struct Design { std::vector factor_names; std::vector sample_names; @@ -69,7 +63,7 @@ operator<<(std::ostream &os, const Design &design); struct mcounts { std::uint32_t n_reads{}; - std::uint32_t n_meth{}; + double n_meth{}; }; [[nodiscard]] inline std::istream & @@ -85,6 +79,18 @@ struct SiteProp { parse(const std::string &line); }; +struct vars_cache { + double p{}; + double a{}; + double b{}; + double lgamma_a{}; + double lgamma_b{}; + double lgamma_a_b{}; + double digamma_a{}; + double digamma_b{}; + double digamma_a_b{}; +}; + struct Regression { static double tolerance; // 1e-3; static double stepsize; // 0.001; @@ -94,11 +100,7 @@ struct Regression { SiteProp props; double max_loglik{}; - // scratch space - std::vector cumul; - std::vector p_v; - std::vector cache; - std::uint32_t max_r_count{}; + std::vector cache; // scratch space [[nodiscard]] std::size_t n_factors() const { @@ -125,6 +127,11 @@ struct Regression { return design.n_samples(); } + [[nodiscard]] std::string + rowname() const { + return props.rowname; + } + void order_samples(const std::vector &ordered_names) { design.order_samples(ordered_names); diff --git a/src/radmeth/radmeth_optimize.cpp b/src/radmeth/radmeth_optimize.cpp index fcbeb1eb..6c6e2621 100644 --- a/src/radmeth/radmeth_optimize.cpp +++ b/src/radmeth/radmeth_optimize.cpp @@ -17,6 +17,7 @@ #include "radmeth_model.hpp" #include +#include // for gsl_sf_psi (digamma) #include #include @@ -36,132 +37,128 @@ get_p(const std::vector &v, const gsl_vector *params) { return logistic(std::inner_product(a, a + std::size(v), params->data, 0.0)); } -static inline auto -set_max_r_count(Regression ®) { - const auto &cumul = reg.cumul; - const auto max_itr = std::max_element( - std::cbegin(cumul), std::cend(cumul), [](const auto &a, const auto &b) { - return std::size(a.r_counts) < std::size(b.r_counts); - }); - reg.max_r_count = std::size(max_itr->r_counts); - // ADS: avoid the realloc that can happen even for resize(smaller_size) - if (reg.max_r_count > std::size(reg.cache)) - reg.cache.resize(reg.max_r_count); -} - -static inline auto -cache_log1p_factors(Regression ®, const double phi) { - const std::size_t max_k = reg.max_r_count; - auto &cache = reg.cache; - for (std::size_t k = 0; k < max_k; ++k) - cache[k] = std::log1p(phi * (k - 1.0)); -} - -static inline auto -cache_dispersion_effect(Regression ®, const double phi) { - const std::size_t max_k = reg.max_r_count; - auto &cache = reg.cache; - for (std::size_t k = 0; k < max_k; ++k) - cache[k] = (k - 1.0) / (1.0 + phi * (k - 1.0)); +static void +get_cache_lgamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.lgamma_a = std::lgamma(a); + cache.lgamma_b = std::lgamma(b); + cache.lgamma_a_b = std::lgamma(a + b); } [[nodiscard]] static double log_likelihood(const gsl_vector *params, Regression ®) { - const auto phi = logistic(gsl_vector_get(params, reg.design.n_factors())); - const auto one_minus_phi = 1.0 - phi; + const auto phi = 1.0 / std::exp(gsl_vector_get(params, reg.n_factors())); const auto n_groups = reg.n_groups(); const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - // ADS: precompute the log1p(phi * (k - 1.0)) values, which are reused for - // each group. - cache_log1p_factors(reg, phi); - const auto &log1p_fact_v = reg.cache; - - double log_lik = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = get_p(groups[g_idx], params); - const auto one_minus_p = 1.0 - p; - - const auto term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - for (std::size_t k = 0; k < std::size(cumul_y); ++k) - log_lik += cumul_y[k] * std::log(term1 + phi * k); - - const auto term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - for (std::size_t k = 0; k < std::size(cumul_d); ++k) - log_lik += cumul_d[k] * std::log(term2 + phi * k); - - const auto &cumul_n = cumul[g_idx].r_counts; - for (std::size_t k = 0; k < std::size(cumul_n); ++k) - log_lik -= cumul_n[k] * log1p_fact_v[k]; + auto &cache = reg.cache; + for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) + get_cache_lgamma(groups[g_idx], params, phi, cache[g_idx]); + + const auto &group_id = reg.design.group_id; + const auto &mc = reg.props.mc; + double ll = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto a = c.a; + const auto b = c.b; + + // clang-format off + ll += ((std::lgamma(y + a) - c.lgamma_a) + + (std::lgamma(n - y + b) - c.lgamma_b) + + (c.lgamma_a_b - std::lgamma(n + a + b))); + // clang-format on } - return log_lik; + + return ll; +} + +[[nodiscard]] static inline double +digamma(const double x) { + return gsl_sf_psi(x); +} + +static void +get_cache_digamma(const std::vector &group, + const gsl_vector *params, const double phi, + vars_cache &cache) { + const auto p = get_p(group, params); + const auto a = p * phi; + const auto b = (1.0 - p) * phi; + cache.p = p; + cache.a = a; + cache.b = b; + cache.digamma_a = digamma(a); + cache.digamma_b = digamma(b); + cache.digamma_a_b = digamma(a + b); } static void gradient(const gsl_vector *params, Regression ®, gsl_vector *output) { - const auto n_factors = reg.design.n_factors(); - const auto phi = logistic(gsl_vector_get(params, n_factors)); - const auto one_minus_phi = 1.0 - phi; + const auto n_factors = reg.n_factors(); + const auto phi = 1.0 / std::exp(gsl_vector_get(params, n_factors)); + + const auto &mc = reg.props.mc; + const auto &matrix = reg.design.matrix; const auto n_groups = reg.n_groups(); const auto &groups = reg.design.groups; - const auto &cumul = reg.cumul; - - auto &p_v = reg.p_v; // ADS: reusing scratch space + auto &cache = reg.cache; // ADS: reusing scratch space for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) - p_v[g_idx] = get_p(groups[g_idx], params); - - cache_dispersion_effect(reg, phi); - const auto &dispersion_effect = reg.cache; // (k-1)/(1 + phi(k-1)) + get_cache_digamma(groups[g_idx], params, phi, cache[g_idx]); - // init output to zero for all factors - gsl_vector_set_all(output, 0.0); + gsl_vector_set_all(output, 0.0); // init output to zero for all factors auto &data = output->data; - double disp_deriv = 0.0; - for (std::size_t g_idx = 0; g_idx < n_groups; ++g_idx) { - const auto p = p_v[g_idx]; - const auto one_minus_p = 1.0 - p; - - double deriv = 0.0; - - const auto denom_term1 = one_minus_phi * p; - const auto &cumul_y = cumul[g_idx].m_counts; - const auto y_lim = std::size(cumul_y); - for (auto k = 0u; k < y_lim; ++k) { - const auto common_factor = cumul_y[k] / (denom_term1 + phi * k); - deriv += common_factor; - disp_deriv += (k - p) * common_factor; - } - - const auto denom_term2 = one_minus_phi * one_minus_p; - const auto &cumul_d = cumul[g_idx].d_counts; - const auto d_lim = std::size(cumul_d); - for (auto k = 0u; k < d_lim; ++k) { - const auto common_factor = cumul_d[k] / (denom_term2 + phi * k); - deriv -= common_factor; - disp_deriv += (k - one_minus_p) * common_factor; - } - - const auto &cumul_n = cumul[g_idx].r_counts; - const auto n_lim = std::size(cumul_n); - for (auto k = 0u; k < n_lim; ++k) - disp_deriv -= cumul_n[k] * dispersion_effect[k]; - - const auto &g = groups[g_idx]; - const auto denom_term1_one_minus_p = denom_term1 * one_minus_p; - for (auto fact_idx = 0u; fact_idx < n_factors; ++fact_idx) { - const auto level = g[fact_idx]; - if (level == 0) - continue; - data[fact_idx] += deriv * (denom_term1_one_minus_p * level); - } + const auto &group_id = reg.design.group_id; + double grad_phi = 0.0; + const auto n_samples = reg.n_samples(); + for (auto i = 0u; i < n_samples; ++i) { + const auto y = mc[i].n_meth; + const auto n = mc[i].n_reads; + + const auto &c = cache[group_id[i]]; + const auto p = c.p; + const auto a = c.a; + const auto b = c.b; + + const auto digamma_a = c.digamma_a; + const auto digamma_b = c.digamma_b; + const auto digamma_y_a = digamma(y + a); + const auto digamma_n_y_b = digamma(n - y + b); + + // grad wrt p + const auto dlogl_dp = + phi * (digamma_y_a - digamma_a - digamma_n_y_b + digamma_b); + const auto dp_delta = p * (1.0 - p); // chain rule: d p / d param + const auto dlogl_delta = dlogl_dp * dp_delta; // grad wrt params + + auto matrix_itr = std::cbegin(matrix[i]); + const auto data_end = data + n_factors; + for (auto data_itr = data; data_itr != data_end; ++data_itr) + *data_itr += (*matrix_itr++) * dlogl_delta; + + const auto digamma_delta = c.digamma_a_b - digamma(n + a + b); + const auto dphi_term1 = p * (digamma_y_a - digamma_a + digamma_delta); + const auto dphi_term2 = + (1.0 - p) * (digamma_n_y_b - digamma_b + digamma_delta); + + grad_phi += dphi_term1 + dphi_term2; } - gsl_vector_set(output, n_factors, disp_deriv * (phi * one_minus_phi)); + + const auto grad_theta = -grad_phi * phi; + gsl_vector_set(output, n_factors, grad_theta); } [[nodiscard]] static double @@ -184,63 +181,14 @@ neg_loglik_and_grad(const gsl_vector *params, void *object, double *loglik_val, neg_gradient(params, object, d_loglik_val); } -static void -get_cumulative(const std::vector &group_id, - const std::uint32_t n_groups, const std::vector &mc, - std::vector &cumul) { - const auto n_cols = std::size(mc); - cumul.clear(); - cumul.resize(n_groups); - - const auto comp_cumul = [&](auto get_value, auto get_vector) { - // phase 1: determine max value for each group - for (auto g_idx = 0u; g_idx < n_groups; ++g_idx) { - std::uint32_t max_v{}; - for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { - if (group_id[c_idx] == g_idx) { - const auto val = get_value(mc[c_idx]); - if (val > max_v) - max_v = val; - } - } - get_vector(cumul[g_idx]).resize(max_v, 0); - } - - // phase 2: fill cumulative counts - for (auto c_idx = 0u; c_idx < n_cols; ++c_idx) { - const auto g_idx = group_id[c_idx]; - const auto val = get_value(mc[c_idx]); - auto &vec = get_vector(cumul[g_idx]); - for (auto i = 0u; i < val; ++i) - ++vec[i]; - } - }; - // call the lambda 3 times for m_counts, r_counts, d_counts - comp_cumul( - [](const mcounts &m) { return m.n_meth; }, - [](cumul_counts &c) -> std::vector & { return c.m_counts; }); - - comp_cumul( - [](const mcounts &m) { return m.n_reads; }, - [](cumul_counts &c) -> std::vector & { return c.r_counts; }); - - comp_cumul( - [](const mcounts &m) { return m.n_reads - m.n_meth; }, - [](cumul_counts &c) -> std::vector & { return c.d_counts; }); -} - void fit_regression_model(Regression &r, std::vector &p_estimates, double &dispersion_estimate) { static constexpr auto init_dispersion_param = -2.5; const auto stepsize = Regression::stepsize; const auto max_iter = Regression::max_iter; - const auto n_groups = r.n_groups(); - get_cumulative(r.design.group_id, n_groups, r.props.mc, r.cumul); - set_max_r_count(r); - - r.p_v.resize(n_groups); + r.cache.resize(n_groups); // make sure scratch space is allocated const std::size_t n_params = r.n_params(); const auto tol = std::sqrt(n_params) * r.n_samples() * Regression::tolerance;