Draft of functions for quantifying power of families for rare variant discovery

.github/ISSUE_TEMPLATE/bug_report.md

@@ -0,0 +1,34 @@
+---
+name: Bug Report
+about: Create a report about something not working as expected.
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+## **Describe the bug**
+A clear and concise description of what the bug is.
+
+## **To reproduce the bug**
+Steps to reproduce the observed issue:
+1. Go to '...'
+2. Click on '....'
+3. Scroll down to '....'
+4. See error
+
+## **Expected behaviour**
+A clear and concise description of what you expected to happen.
+
+## **Screenshots and data**
+If applicable, add input/output data (files or cloud location) or screenshots to help explain your problem.
+
+## **Additional context**
+Add any other context about the problem here.
+
+## **Proposed troubleshooting steps**
+Describe how you're going to try to fix it!
+- [ ] This can be a task list
+
+## **Help requests**
+Describe any specific aspects of the problem/solution where help from a teammate is required. Tag them so they're aware.

.github/ISSUE_TEMPLATE/feature-addition.md

@@ -0,0 +1,26 @@
+---
+name: Feature Addition
+about: Suggest an idea for this project.
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+## **Is your feature request related to a problem? Please describe.**
+A clear and concise description of what the problem is. Ex. Currently we cannot test ___, the way ___ is done is disorganized [...]
+
+## **Describe the solution you'd like**
+A clear and concise description of what you want to happen, in plain english.
+
+## **Describe the plan**
+A clear and concise description of the path forward as you see it. Describe any uncertainties or alternative solutions that need to be compared.
+
+### **TODO**
+- [ ] Break the plan down into specifics
+
+## **Additional context**
+Add any other context or screenshots about the feature request here.
+
+## **Help requests**
+Describe any specific aspects of the problem/solution where help from a teammate is required. Tag them so they're aware.

.github/pull_request_template.md

@@ -0,0 +1,25 @@
+# Overview of Pull Request
+A brief description of _why_ this change is necessary (couple of sentences/short paragraph).
+
+## Release notes
+A brief description of what this change is in plain english. What does this change add, fix, or address? (couple of sentences/short paragraph)
+
+## Closes the following issues
+ -  _[Links to any relevant GitHub Issues]_
+
+## Pre-review checklist
+Make sure all these boxes are checked before tagging a reviewer!
+
+- [ ] The PR title is a concise, present-tense summary of the change
+- [ ] The PR is linked to any relevant GitHub Issues, if they exist
+- [ ] The PR is against the intended branch, and is mergeable
+- [ ] The code is sufficiently tested (unit tests and real world experiments, where applicable)
+- [ ] The code is linted and meets style guidelines
+- [ ] The changelog has been updated
+- [ ] I reviewed the PR myself before requesting a review from others
+
+## Pre-release notes
+* External dependencies affected
+    - [ ] _List any dependencies that may be affected by this change or need to be updated to align with this change_
+* Post-release tasks
+    - [ ] _List any rake tasks or other TODOs that need to be sorted after this change is released_

.github/workflows/main.yaml

@@ -0,0 +1,37 @@
+name: ContiniousIntegrationChecks
+run-name: ${{ github.actor }} triggered CI run with event ${{ github.event_name }} for branch ${{ github.ref }}
+on:
+  push:
+    branches: [ main, ci_add]
+  pull_request:
+    branches_ignore: []
+permissions:
+  id-token: write
+  contents: read # This is required for actions/checkout
+jobs:
+  fullCI:
+    defaults:
+      run:
+        shell: bash -l {0}
+    name: CI test
+    runs-on: ubuntu-latest
+    steps:
+      - run: echo "The job was automatically triggered by a ${{ github.event_name }} event."
+      - run: echo "The name of your branch is ${{ github.ref }} and your repository is ${{ github.repository }}."
+      - run: echo "The job's status is ${{ job.status }}."
+      - name: Git clone the repository
+        uses: actions/checkout@v4
+      - name: Set up R env
+        uses: r-lib/actions/setup-r@v2
+      - name: Set up conda
+        uses: conda-incubator/setup-miniconda@v3
+        with:
+          auto-activate-base: false
+          python-version: 3.12
+          channels: bioconda, conda-forge, defaults
+      - name: Conda create env and install dependencies
+        shell: bash -l {0}
+        run: |
+          conda create -n test_r r-base r-devtools r-testthat
+          conda activate test_r
+          Rscript -e "testthat::test_local()" 

DESCRIPTION

@@ -1,22 +1,23 @@
-Package: seqbio.variant.scoring
-Title: Score evidence of rare variants' association with disease status in a family, based on relationships of individuals
+Package: KinformR
+Title: Relationship-informed pedigree and variant scoring
 Version: 0.1.0
 Authors@R: 
     person("Cameron", "Nugent", , "cam.nugent@sequencebio.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-1135-2605"))
 Description: 
-    This R package is designed to score rare variants, assigning values based on 
-    the disease status of individuals, the presence or absence of a rare variant 
-    in those individuals, and their pairwise coefficients of relatedness. 
-    The package uses a custom formula to assign value to a variant that gives 
-    more weight to shared variants common to distantly related affected 
-    individuals. The variant status for unaffected individuals can optionally 
-    be considered as well, with the highest scoring values being given to 
-    closely related individuals that do not share a variant of interst. 
-    Since variants can be incompletely penetrant, the scoring can be based 
-    solely on the affected individuals, or the weight of unaffected evidence 
-    can be customized.
+    The KinformR R package is meant to aid in comparative evaluation of families 
+    and candidate variants in rare-variant association studies. The package can be used for 
+    two methodologically overlapping but distinct purposes. First, the prior to any genetic or genomic 
+    evaluation, evaluation of relative detection power of pedigrees, can direct recruitment 
+    efforts by showing which unsampled individuals would be the most meaningful additions to a study. 
+    Second, after sequencing and analysis,  variants based on association with disease status 
+    and familial relationships of individuals, aids in variant prioritization.
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.2
+VignetteBuilder: knitr
+Suggests:
+    devtools,
+    testthat,
+    knitr

LICENSE

@@ -1,2 +1,21 @@
-YEAR: 2025
-COPYRIGHT HOLDER: seqbio.variant.scoring authors
+# MIT License
+
+Copyright (c) 2025 KinformR Sequence Bioinformatics Inc.
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

NAMESPACE

@@ -1,12 +1,17 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method(sum,fam.scores)
+export(add.fam.scores)
 export(assign.status)
 export(build.relation.dict)
 export(calc.rv.score)
 export(encode.rows)
+export(ibd)
+export(penetrance)
 export(read.indiv)
+export(read.pedigree)
 export(read.relation.mat)
 export(read.var.table)
+export(score)
 export(score.fam)
+export(score.pedigree)
 export(score.variant.status)

R/encoding.R

@@ -1,12 +1,12 @@
 
 
 #' Take a disease status and a genetic variant and determine which category the combo falls in.
-#' A_c = Affected individual with ALT variant
-#' A_i = Affected individual without ALT variant
-#' U_c = Unaffected individual without ALT variant
-#' U_i = Unaffected individual with ALT variant
+#' A.c = Affected individual with ALT variant
+#' A.i = Affected individual without ALT variant
+#' U.c = Unaffected individual without ALT variant
+#' U.i = Unaffected individual with ALT variant
 #' If theoretical.max = TRUE the true variant statuses are ignored and all
-#' affected/unaffected are assigned A_c and U_c respectively.
+#' affected/unaffected are assigned A.c and U.c respectively.
 #' These encoding can then be used show what a family's max score would be.
 #'
 #' @param status Disease status of an individual. A = affected, U = unaffected.
@@ -16,41 +16,44 @@
 #' Default is FALSE.
 #' @return a string
 #' @examples
-#' assign.status("A", "0/1") == "A_c"
-#' assign.status("A", "0|0") == "A_i"
-#' assign.status("U", 1) == "U_i"
-#' assign.status("U", "0|0") =="U_c"
+#' assign.status("A", "0/1") == "A.c"
+#' assign.status("A", "0|0") == "A.i"
+#' assign.status("U", 1) == "U.i"
+#' assign.status("U", "0|0") =="U.c"
 #' @export
 assign.status <- function(status, variant,  theoretical.max=FALSE){
+
+  var.err<-"Incompatible variant value! Supported encodings are: '0' '1' '0/0' '0/1' '0|0' '0|1'"
   if(status == "A"){
     if(theoretical.max){
-      return("A_c")
+      return("A.c")
     }
-    else if(variant == "0/1" || variant == "1/1" || variant == "1" || variant == "0|1" || variant == "1|0"|| variant == "1|1"  ){
-      return("A_c")
+    #NOTE - Once in a while 1/0 genotypes crop up; also 0/2 etc. if derived from multi-allelics. This edge case not covered at present.
+    else if(variant == "0/1" || variant == "1/1" || variant == "1" || variant == "0|1" || variant == "1|0" || variant == "1|1"  ){
+      return("A.c")
     }else if (variant == "0/0" || variant == "0" || variant == "0|0" ){
-      return("A_i")
+      return("A.i")
     }else{
-      return("unk")
+      stop(var.err)
     }
   }else if (status == "U"){
     if(theoretical.max){
-      return("U_c")
-    } else if(variant == "0/1" || variant == "1/1" || variant == "1" || variant == "0|1" || variant == "1|0"|| variant == "1|1"  ){
-      return("U_i")
+      return("U.c")
+    } else if(variant == "0/1" || variant == "1/1" || variant == "1" || variant == "0|1" || variant == "1|0" || variant == "1|1"  ){
+      return("U.i")
     }else if (variant == "0/0" || variant == "0" || variant == "0|0" ){
-      return("U_c")
+      return("U.c")
     }else{
-      return("unk")
+      stop(var.err)
     }
   }else{
-    return("unk")
+    stop("Status must be one of: U or A")
   }
 }
 
 #' Take the dataframe with variants and status and determine which indivudals
 #' are scored correctly and which are scored incorrectly.
-#' Assign an A_c, A_i, U_c, U_i, unk
+#' Assign an A.c, A.i, U.c, U.i, unk
 #'
 #' Variants can be encoded as binary (0 or 1, genotypes 0/0 or 0/1, or phased genotypes 0|0 0|1).
 #' Note the program assumes alt is the disease allele. homozygous alts are allowed.
@@ -61,7 +64,7 @@ assign.status <- function(status, variant,  theoretical.max=FALSE){
 #' TODO - switch to numbers 1-4 and -1?
 #' @param indiv.df A dataframe with the format:
 #' name	         status	variant
-#' MS-4107-1001      A      0/1
+#' MS-5678-1001      A      0/1
 #' @param theoretical.max Should the theoretical maxima be returned instead of the observed values?
 #' When true, the scoring assumes correct variant-status pair for each individual.
 #' Default is FALSE.
@@ -95,15 +98,13 @@ return(indiv.df)
 #' @param drop.unrelated Should unrelated (-1) relationships be dropped? Default = TRUE.
 #'
 #' @return A list with the categorized relationship/variant information.
-#' @examples
-#' #TODO - add
 #' @export
 build.relation.dict <- function( mat.row, name.stat.dict, drop.unrelated=TRUE){
   indiv.rels = list(
-    "A_c" = c(),
-    "A_i" = c(),
-    "U_c" = c(),
-    "U_i" = c()
+    "A.c" = c(),
+    "A.i" = c(),
+    "U.c" = c(),
+    "U.i" = c()
   )
 
   for(i in seq_along(mat.row)){
@@ -124,10 +125,9 @@ build.relation.dict <- function( mat.row, name.stat.dict, drop.unrelated=TRUE){
 #' For each row, generate the encoded data for scoring.
 #' @param relation.mat The relationship matrix for all pairwise combinations of individuals.
 #' @param status.df The ID, status, and genotypes for each individual.
+#' @param ... Additional arguments to be passed between methods.
 #' @return A dictionary with the per-individual relationship lists.
 #' One value for each row of the matrix.
-#' @examples
-#' #TODO - add
 #' @export
 encode.rows <- function(relation.mat, status.df, ...){
 

R/io.R

@@ -4,10 +4,10 @@
 #'
 #' @param fname A file name, expected format of contents is:
 #' name	         status	variant
-#' MS-4107-1001      A      0/1
+#' MS-5678-1001      A      0/1
 #' @return A data frame.
 #' @examples
-#' tsv.name1 <-system.file('extdata/7003_notch3.tsv', package = 'seqbio.variant.scoring')
+#' tsv.name1 <-system.file('extdata/1234_ex2.tsv', package = 'KinformR')
 #' id.df1 <- read.indiv(tsv.name1)
 #' @export
 read.indiv <- function(fname){
@@ -24,7 +24,7 @@ read.indiv <- function(fname){
 #' @param fname The file with the relationship matrix information.
 #' @return A matrix with the relationships and individual ids as rownames and colnames.
 #' @examples
-#' mat.name1 <-system.file('extdata/7003_notch3.mat', package = 'seqbio.variant.scoring')
+#' mat.name1 <-system.file('extdata/1234_ex2.mat', package = 'KinformR')
 #' mat1 <- read.relation.mat(mat.name1)
 #' @export
 read.relation.mat <- function(fname){
@@ -49,15 +49,15 @@ read.relation.mat <- function(fname){
 #'
 #'
 #' @param fname A file name, expected format of contents is:
-#' #CHROM  POS       REF  ALT  MS-4107-1001_A  MS-4107-1002_U  ...
+#' #CHROM  POS       REF  ALT  MS-5678-1001_A  MS-5678-1002_U  ...
 #' chr3    46203838  G    A    0/1             0/0      ...
 #' @return A dataframe.
 #' Data will be worked into a data frame with format.
 #' name	         status	variant
-#' MS-4107-1001      A      0/1
+#' MS-5678-1001      A      0/1
 #' @examples
-#' ex.infile <-system.file('extdata/example_vcf_extract_4107.tsv',
-#'                          package = 'seqbio.variant.scoring')
+#' ex.infile <-system.file('extdata/example_vcf_extract_5678.tsv',
+#'                          package = 'KinformR')
 #' read.var.table(ex.infile)
 #' @export
 read.var.table <- function(fname){