addReadStats/calcReadStats

fmicompbio / footprintR

Tools for working with single-molecule footprinting data

https://fmicompbio.github.io/footprintR/

Other

0 stars 0 forks source link

addReadStats/calcReadStats #17

Closed csoneson closed 3 weeks ago

csoneson commented 1 month ago

[x] add argument to addReadStats to name the slot (defaults to QC): done in 7487ac87
[x] append metadata to SimpleList returned by calcReadStats (in calcReadStats), containing parameters: done in 55f5d079 (metadata(calcReadStats(...))
[x] keep position filtering steps (default turned off) in calcReadStats: also in 55f5d079
[x] make sure calcReadStats works with NaArray: was always working
[x] add quality score: see next item
[x] add read length/alignment length/aligned fraction: added in 2da3f017 (C++) and 4f8d8f8b (R): is returned by readModBam in colData(se)$read_info (SimpleList with one element per sample)

mbstadler commented 1 month ago

Preparing for the above, I had to update read_modbam_cpp (2da3f017):

add calculate_aligned_bases helper function
fix non-persistent qseq bug
add read_length and aligned_length return values
return read-level data as a data.frame (instead of repeated for each modified base)
return read-level data only for reads with at least one modified base
store bamdata->core.l_qseq to avoid repeated pointer calculations
update documentation
adapt unit tests

mbstadler commented 1 month ago

And as a result, I updated readModBam and its unit tests to work with the new return value of read_modbam_cpp (4f8d8f8)

mbstadler commented 1 month ago

I think this is all done now - however only if the data is read using readModBam. We should discuss if information like qscore and aligned_length can also be obtained when using readModkitExtract, or if they are needed in that case.

mbstadler commented 1 month ago

Here is what we are currently extracting in terms of read-level information in readModBam:

modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "footprintR")
se <- readModBam(bamfiles = modbamfile, regions = "chr1:6940000-6955000", modbase = "a", verbose = TRUE)
se$read_info$s1
                                  read_id  qscore read_length aligned_length aligned_fraction
1 s1-233e48a7-f379-4dcf-9270-958231125563 14.1428       20058          14801        0.7379101
2 s1-d52a5f6a-a60a-4f85-913e-eada84bfbfb9 16.0127       11305          11214        0.9919505
3 s1-92e906ae-cddb-4347-a114-bf9137761a8d 20.3082       12277          12227        0.9959273

Here is what information could possibly be extracted by readModkitExtract:

qscore: not possible (only have the qscore for individual modified bases)
read_length: doable
aligned_length: not possible (only have read_length and soft-masked lengths, but lack information about indels)
aligned_fraction: not possible (relies on aligned_length)

mbstadler commented 3 weeks ago

Given the above summary, it seems unlikely that we can make readModkitExtract extract the same information as readModBam. I am closing this issue (feel free to reopen it if necessary).