magnusdv / forrel Goto Github PK

View Code? Open in Web Editor NEW

9.0 9.0 0.0 1.53 MB

Forensic pedigree analysis and relatedness inference

License: GNU General Public License v2.0

R 100.00%

forrel's People

Contributors

Stargazers

Watchers

forrel's Issues

readFam contd

Previous example now runs. Here's an example probably best ignored.
readFam fails for an old version of Familias (3.1.9), but works if the file is saved as most
recent Familias version (3.2.9). Here are the runs in turn:

x = forrel::readFam("http://familias.name/BookKEP/Solution3_5.fam") # version 3.1.9, fails
#> Familias version: 3.1.9 (DVI)
#> Read DVI: Yes
#> 
#> Number of individuals (excluding 'extras'): 100
#>   Individual 'Sister1': Genotypes for 22 markers read
#>   Individual 'Sister2': Genotypes for 22 markers read
#>   Individual 'Sister3': Genotypes for 22 markers read
#>   Individual 'Sister4': Genotypes for 22 markers read
#>   Individual 'Sister5': Genotypes for 22 markers read
#>   Individual 'Sister6': Genotypes for 22 markers read
#>   Individual 'Sister7': Genotypes for 22 markers read
#>   Individual 'Sister8': Genotypes for 22 markers read
#>   Individual 'Sister9': Genotypes for 22 markers read
#>   Individual 'Sister10': Genotypes for 22 markers read
#>   Individual 'Sister11': Genotypes for 22 markers read
#>   Individual 'Sister12': Genotypes for 22 markers read
#>   Individual 'Sister13': Genotypes for 22 markers read
#>   Individual 'Sister14': Genotypes for 22 markers read
#>   Individual 'Sister15': Genotypes for 22 markers read
#>   Individual 'Sister16': Genotypes for 22 markers read
#>   Individual 'Sister17': Genotypes for 22 markers read
#>   Individual 'Sister18': Genotypes for 22 markers read
#>   Individual 'Sister19': Genotypes for 22 markers read
#>   Individual 'Sister20': Genotypes for 22 markers read
#>   Individual 'Sister21': Genotypes for 22 markers read
#>   Individual 'Sister22': Genotypes for 22 markers read
#>   Individual 'Sister23': Genotypes for 22 markers read
#>   Individual 'Sister24': Genotypes for 22 markers read
#>   Individual 'Sister25': Genotypes for 22 markers read
#>   Individual 'Sister26': Genotypes for 22 markers read
#>   Individual 'Sister27': Genotypes for 22 markers read
#>   Individual 'Sister28': Genotypes for 22 markers read
#>   Individual 'Sister29': Genotypes for 22 markers read
#>   Individual 'Sister30': Genotypes for 22 markers read
#>   Individual 'Sister31': Genotypes for 22 markers read
#>   Individual 'Sister32': Genotypes for 22 markers read
#>   Individual 'Sister33': Genotypes for 22 markers read
#>   Individual 'Sister34': Genotypes for 22 markers read
#>   Individual 'Sister35': Genotypes for 22 markers read
#>   Individual 'Sister36': Genotypes for 22 markers read
#>   Individual 'Sister37': Genotypes for 22 markers read
#>   Individual 'Sister38': Genotypes for 22 markers read
#>   Individual 'Sister39': Genotypes for 22 markers read
#>   Individual 'Sister40': Genotypes for 22 markers read
#>   Individual 'Sister41': Genotypes for 22 markers read
#>   Individual 'Sister42': Genotypes for 22 markers read
#>   Individual 'Sister43': Genotypes for 22 markers read
#>   Individual 'Sister44': Genotypes for 22 markers read
#>   Individual 'Sister45': Genotypes for 22 markers read
#>   Individual 'Sister46': Genotypes for 22 markers read
#>   Individual 'Sister47': Genotypes for 22 markers read
#>   Individual 'Sister48': Genotypes for 22 markers read
#>   Individual 'Sister49': Genotypes for 22 markers read
#>   Individual 'Sister50': Genotypes for 22 markers read
#>   Individual 'Sister51': Genotypes for 22 markers read
#>   Individual 'Sister52': Genotypes for 22 markers read
#>   Individual 'Sister53': Genotypes for 22 markers read
#>   Individual 'Sister54': Genotypes for 22 markers read
#>   Individual 'Sister55': Genotypes for 22 markers read
#>   Individual 'Sister56': Genotypes for 22 markers read
#>   Individual 'Sister57': Genotypes for 22 markers read
#>   Individual 'Sister58': Genotypes for 22 markers read
#>   Individual 'Sister59': Genotypes for 22 markers read
#>   Individual 'Sister60': Genotypes for 22 markers read
#>   Individual 'Sister61': Genotypes for 22 markers read
#>   Individual 'Sister62': Genotypes for 22 markers read
#>   Individual 'Sister63': Genotypes for 22 markers read
#>   Individual 'Sister64': Genotypes for 22 markers read
#>   Individual 'Sister65': Genotypes for 22 markers read
#>   Individual 'Sister66': Genotypes for 22 markers read
#>   Individual 'Sister67': Genotypes for 22 markers read
#>   Individual 'Sister68': Genotypes for 22 markers read
#>   Individual 'Sister69': Genotypes for 22 markers read
#>   Individual 'Sister70': Genotypes for 22 markers read
#>   Individual 'Sister71': Genotypes for 22 markers read
#>   Individual 'Sister72': Genotypes for 22 markers read
#>   Individual 'Sister73': Genotypes for 22 markers read
#>   Individual 'Sister74': Genotypes for 22 markers read
#>   Individual 'Sister75': Genotypes for 22 markers read
#>   Individual 'Sister76': Genotypes for 22 markers read
#>   Individual 'Sister77': Genotypes for 22 markers read
#>   Individual 'Sister78': Genotypes for 22 markers read
#>   Individual 'Sister79': Genotypes for 22 markers read
#>   Individual 'Sister80': Genotypes for 22 markers read
#>   Individual 'Sister81': Genotypes for 22 markers read
#>   Individual 'Sister82': Genotypes for 22 markers read
#>   Individual 'Sister83': Genotypes for 22 markers read
#>   Individual 'Sister84': Genotypes for 22 markers read
#>   Individual 'Sister85': Genotypes for 22 markers read
#>   Individual 'Sister86': Genotypes for 22 markers read
#>   Individual 'Sister87': Genotypes for 22 markers read
#>   Individual 'Sister88': Genotypes for 22 markers read
#>   Individual 'Sister89': Genotypes for 22 markers read
#>   Individual 'Sister90': Genotypes for 22 markers read
#>   Individual 'Sister91': Genotypes for 22 markers read
#>   Individual 'Sister92': Genotypes for 22 markers read
#>   Individual 'Sister93': Genotypes for 22 markers read
#>   Individual 'Sister94': Genotypes for 22 markers read
#>   Individual 'Sister95': Genotypes for 22 markers read
#>   Individual 'Sister96': Genotypes for 22 markers read
#>   Individual 'Sister97': Genotypes for 22 markers read
#>   Individual 'Sister98': Genotypes for 22 markers read
#>   Individual 'Sister99': Genotypes for 22 markers read
#>   Individual 'Sister100': Genotypes for 22 markers read
#> 
#> Number of pedigrees: 0
#> 
#> Database: Norwegian database
#> Number of loci: 22
#>   D3S1358: 11 alleles, model = equal, rate = 0 (unisex)
#>   TH01: 8 alleles, model = equal, rate = 0 (unisex)
#>   D21S11: 19 alleles, model = equal, rate = 0 (unisex)
#>   D18S51: 20 alleles, model = equal, rate = 0 (unisex)
#>   PENTA_E: 18 alleles, model = equal, rate = 0 (unisex)
#>   D5S818: 9 alleles, model = equal, rate = 0 (unisex)
#>   D13S317: 9 alleles, model = equal, rate = 0 (unisex)
#>   D7S820: 15 alleles, model = equal, rate = 0 (unisex)
#>   D16S539: 9 alleles, model = equal, rate = 0 (unisex)
#>   CSF1PO: 10 alleles, model = equal, rate = 0 (unisex)
#>   PENTA_D: 16 alleles, model = equal, rate = 0 (unisex)
#>   VWA: 12 alleles, model = equal, rate = 0 (unisex)
#>   D8S1179: 10 alleles, model = equal, rate = 0 (unisex)
#>   TPOX: 8 alleles, model = equal, rate = 0 (unisex)
#>   FGA: 21 alleles, model = equal, rate = 0 (unisex)
#>   D12S391: 19 alleles, model = equal, rate = 0 (unisex)
#>   SE33: 50 alleles, model = equal, rate = 0 (unisex)
#>   D10S1248: 7 alleles, model = equal, rate = 0 (unisex)
#>   D1S1656: 14 alleles, model = equal, rate = 0 (unisex)
#>   D2S1338: 11 alleles, model = equal, rate = 0 (unisex)
#>   D22S1045: 8 alleles, model = equal, rate = 0 (unisex)
#>   D2S441: 9 alleles, model = equal, rate = 0 (unisex)
#> 
#> *** Reading DVI section ***
#> Unidentified persons: 100
#> Error in asFamiliasPedigree(as.character(id), 0, 0, as.integer(sex)): 'list' object cannot be coerced to type 'integer'

^{Created on 2020-02-12 by the reprex package (v0.3.0)}

x =forrel::readFam("http://familias.name/BookKEP/Solution3_5-3.2.9.fam") 
#> Familias version: 3.2.9 (DVI)
#> Read DVI: Yes
#> 
#> Number of individuals (excluding 'extras'): 0
#> 
#> Number of pedigrees: 0
#> 
#> Database: Norwegian database
#> Number of loci: 22
#>   D3S1358: 11 alleles, model = equal, rate = 0 (unisex)
#>   TH01: 8 alleles, model = equal, rate = 0 (unisex)
#>   D21S11: 19 alleles, model = equal, rate = 0 (unisex)
#>   D18S51: 20 alleles, model = equal, rate = 0 (unisex)
#>   PENTA_E: 18 alleles, model = equal, rate = 0 (unisex)
#>   D5S818: 9 alleles, model = equal, rate = 0 (unisex)
#>   D13S317: 9 alleles, model = equal, rate = 0 (unisex)
#>   D7S820: 15 alleles, model = equal, rate = 0 (unisex)
#>   D16S539: 9 alleles, model = equal, rate = 0 (unisex)
#>   CSF1PO: 10 alleles, model = equal, rate = 0 (unisex)
#>   PENTA_D: 16 alleles, model = equal, rate = 0 (unisex)
#>   VWA: 12 alleles, model = equal, rate = 0 (unisex)
#>   D8S1179: 10 alleles, model = equal, rate = 0 (unisex)
#>   TPOX: 8 alleles, model = equal, rate = 0 (unisex)
#>   FGA: 21 alleles, model = equal, rate = 0 (unisex)
#>   D12S391: 19 alleles, model = equal, rate = 0 (unisex)
#>   SE33: 50 alleles, model = equal, rate = 0 (unisex)
#>   D10S1248: 7 alleles, model = equal, rate = 0 (unisex)
#>   D1S1656: 14 alleles, model = equal, rate = 0 (unisex)
#>   D2S1338: 11 alleles, model = equal, rate = 0 (unisex)
#>   D22S1045: 8 alleles, model = equal, rate = 0 (unisex)
#>   D2S441: 9 alleles, model = equal, rate = 0 (unisex)
#> 
#> *** Reading DVI section ***
#> Unidentified persons: 100
#>   Sister1
#>   Sister2
#>   Sister3
#>   Sister4
#>   Sister5
#>   Sister6
#>   Sister7
#>   Sister8
#>   Sister9
#>   Sister10
#>   Sister11
#>   Sister12
#>   Sister13
#>   Sister14
#>   Sister15
#>   Sister16
#>   Sister17
#>   Sister18
#>   Sister19
#>   Sister20
#>   Sister21
#>   Sister22
#>   Sister23
#>   Sister24
#>   Sister25
#>   Sister26
#>   Sister27
#>   Sister28
#>   Sister29
#>   Sister30
#>   Sister31
#>   Sister32
#>   Sister33
#>   Sister34
#>   Sister35
#>   Sister36
#>   Sister37
#>   Sister38
#>   Sister39
#>   Sister40
#>   Sister41
#>   Sister42
#>   Sister43
#>   Sister44
#>   Sister45
#>   Sister46
#>   Sister47
#>   Sister48
#>   Sister49
#>   Sister50
#>   Sister51
#>   Sister52
#>   Sister53
#>   Sister54
#>   Sister55
#>   Sister56
#>   Sister57
#>   Sister58
#>   Sister59
#>   Sister60
#>   Sister61
#>   Sister62
#>   Sister63
#>   Sister64
#>   Sister65
#>   Sister66
#>   Sister67
#>   Sister68
#>   Sister69
#>   Sister70
#>   Sister71
#>   Sister72
#>   Sister73
#>   Sister74
#>   Sister75
#>   Sister76
#>   Sister77
#>   Sister78
#>   Sister79
#>   Sister80
#>   Sister81
#>   Sister82
#>   Sister83
#>   Sister84
#>   Sister85
#>   Sister86
#>   Sister87
#>   Sister88
#>   Sister89
#>   Sister90
#>   Sister91
#>   Sister92
#>   Sister93
#>   Sister94
#>   Sister95
#>   Sister96
#>   Sister97
#>   Sister98
#>   Sister99
#>   Sister100
#> 
#> Reference families: 100
#>   Family1 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family2 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family3 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family4 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family5 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family6 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family7 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family8 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family9 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family10 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family11 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family12 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family13 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family14 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family15 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family16 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family17 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family18 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family19 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family20 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family21 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family22 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family23 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family24 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family25 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family26 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family27 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family28 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family29 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family30 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family31 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family32 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family33 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family34 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family35 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family36 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family37 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family38 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family39 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family40 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family41 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family42 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family43 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family44 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family45 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family46 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family47 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family48 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family49 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family50 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family51 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family52 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family53 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family54 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family55 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family56 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family57 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family58 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family59 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family60 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family61 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family62 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family63 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family64 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family65 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family66 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family67 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family68 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family69 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family70 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family71 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family72 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family73 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family74 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family75 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family76 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family77 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family78 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family79 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family80 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family81 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family82 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family83 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family84 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family85 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family86 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family87 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family88 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family89 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family90 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family91 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family92 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family93 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family94 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family95 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family96 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family97 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family98 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family99 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#>   Family100 (3 persons, 2 pedigrees)
#>     Reference pedigree
#>     Missing person
#> 
#> Converting to `ped` format
#> *** Finished DVI section ***

^{Created on 2020-02-12 by the reprex package (v0.3.0)}

Saving and restoring configuration / workspace

Say the user wants to save the configuration of the application so that the next time they want to run the calculations they don't need to import everything again. This would be beneficial for two reasons:

It eliminates repetition and thus errors.
It allows us to add more advanced and extensive options to the program without fear of overwhelming users with settings.

Right now, shiny v0.14 supports bookmarking. Bookmarking makes it very easy for us to save user input (it's really just calling two functions). However, the mechanism that shiny implements for saving is not conventional for desktop applications. Basically, bookmarking works like this:

The user does some work and configures the program to their liking.
The user hits the Bookmark button and the shiny server saves the application state to disk. Then it provides the user with a link.
The user must save the link somewhere to regain access to the configuration they left behin.

Right now, bookmarking has a few shortcomings:

It's not as straightforward to save input files as it is to save simple inputs like checkboxes.
- This could be solved by implementing custom bookmark/restore rules and/or serialisation protocols. See [1].
The application is not fully reactive, meaning that the output at a given time cannot be derived from the input at that moment. In our case, this is due to the exclusion power being calculated after the user hits the Calculate exclusion power button instead of calculating it on the fly when input changes.
- This could also be solved through [1].
Saving multiple times will result in multiple R datastores being saved to disk. This could end up consuming a lot of space.
- See note at the end.
Because R datastores are saved to the inst folder, saving will fail if the user has installed the package in a folder that requires admin access for writing. This might be the case when forrel is used in corporate or lab environments.
For the same reason, updating forrel will likely result in total loss of saved workspaces as the inst folder will be replaced.
- This can be probably solved in a shady way by providing an option to download a .zip file of all saved workspaces and then having users restore it when they upgrade to a new version.
Restoring a workspace requires pasting a URL in the browser. This in turn presents two problems:
1. RStudio users will not typically open their apps in the browser but rather in a window inside RStudio that does not allow them to input a URL manually.
  - Maybe there is an API that allows restoring an arbitrary datastore from its hash. This needs to be researched.
2. The shiny server port might change without notice rendering all URLs unusable and thus saved data inaccessible.
  - See note at the end.
Because shiny runs on a browser and files are uploaded rather than accessed through the filesystem, if users make changes to database, pedigree or reference files, these will not be reflected in the GUI project or workspace if the user does not re-upload them manually.
- There is no straightforward way to fix this. Browser sandboxing prevents random access to the filesystem. Maybe the latests set of web standards has defined an API for this but in any case it would require major work to work with shiny. Anyway, maybe this is not so bad as it seems to me -- users may not be as confused as I think they'll be if their data is not updated when changing the source files.

Other forms of saving application state are not straightforward either. Some of the prior problems still apply, namely those that require interaction with the filesystem. Also, because of the reactive programming paradigm of shiny apps, saving/restoring state is not trivial for some input objects.

Note at the end: one way to solve these two problems at once can be to not show the user the URL after saving (can be done by overriding onBookmarked). Instead, when users want to restore a workspace, they would be presented with a modal containing a list of all saved workspaces. This list can be obtained by listing the contents of the directory in which datastores are saved (inst/exclusionPowerUI/shiny_store). Furthermore, we could allow for only one workspace to be saved at a time, thereby also eliminating the problem with too many saved workspaces filling up space.

References

[1] http://shiny.rstudio.com/articles/advanced-bookmarking.html

`profileSim` for list of singletons

I don't understand why I always get the same alleles for A and B , as below

library(forrel)
#> Loading required package: pedtools
a = singleton("A")
m = marker(a, alleles = 1:1000)
a = setMarkers(a, m)
b = singleton("B")
m = marker(b, alleles = 1:1000)
b = setMarkers(b, m)
U = list(a, b)
profileSim(U,  N = 1,  ids = c("A", "B"), seed = 1234)
#> Preparing parallelisation using 3 cores
#> Preparing parallelisation using 3 cores
#> [[1]]
#> [[1]][[1]]
#>  id fid mid sex     <1>
#>   A   *   *   1 495/860
#> 
#> [[1]][[2]]
#>  id fid mid sex     <1>
#>   B   *   *   1 495/860

^{Created on 2020-08-23 by the reprex package (v0.3.0)}

`exclusionPower` and `founderInbreeding`

It appears that founder inbreeding must be specified for the true alternative in
exclusionPower as the below example shows. Perhaps a clarification is in place
in the documentation

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(nch = 1, sex = 2)
true = list(singleton(id = 1), singleton(id = 3, sex = 2))
m1 = marker(claim)
claim = setMarkers(claim, m1)
f = 1/16
founderInbreeding(claim, "1") = f
# founderInbreeding(true[[1]], "1") = f # gives correct result.
exclusionPower(claim, true, ids = c(1,3), plot = F, verbose = F)
#> Total EP: 0.125 
#> Markers inconsistent with pedigree: 0  
#> Markers with potential mismatch: 1 (1) 
#> Expected number of mismatches: 0.125
# exclusionPower gives the same result, ep = 0.125, for f < 1
# and ep = 0.625 for f = 1. I get
p = 0.5
(f*p+(1-f)*p^2)*(1-p)^2 + (f*(1-p)+(1-f)*(1-p)^2)*p^2
#> [1] 0.1328125

^{Created on 2020-06-15 by the reprex package (v0.3.0)}

`markers` problem in `missingPersonIP`

Great functions: missingPersonIP and missingPersonEP.
So far only one problem detected:

# Four siblings; the fourth is missing
library(forrel, quietly = T)
x = nuclearPed(4)

# Remaining sibs typed with 5 triallelic markers
x = markerSim(x, N = 5, ids = 3:5, alleles = 1:3, seed = 123, verbose = FALSE)

# Compute exclusion power statistics
missingPersonIP(x, missing = 6, markers = 1,
                nsim = 5, threshold = c(10, 100))
#> Simulating 5 profiles...done
#> Computing likelihood ratios...
#> Error in `rownames<-`(`*tmp*`, value = markers): attempt to set 'rownames' on an object with no dimensions

^{Created on 2019-10-26 by the reprex package (v0.3.0)}

`checkPairwise()` fails when exactly two individuals are typed

library(forrel, quietly = T)
x = markerSim(nuclearPed(), N = 10, seed = 123, verbose = F)

# OK with all 3:
checkPairwise(x, plot = F)
#>   id1 id2  N        k0        k1         k2 kappa0 kappa1 kappa2       LR
#> 1   1   2 10 0.9117348 0.0000000 0.08826517      1      0      0 1.088232
#> 2   1   3 10 0.0000000 0.7999559 0.20004411      0      1      0 1.223059
#> 3   2   3 10 0.0000000 0.7999559 0.20004411      0      1      0 1.223059

# Fail:
y = setAlleles(x, ids = 1, alleles = 0)
checkPairwise(y)
#> Error in fix.by(by.y, y): 'by' must match numbers of columns

^{Created on 2022-02-23 by the reprex package (v2.0.1)}

Familias2ped

Here's a problem with Familias2ped:

library(Familias, quietly = TRUE)
persons <- c("Mother", "Child", "Man")
sex <- c("female", "female", "male")
ped2 <- FamiliasPedigree(id = persons, dadid=c(NA,NA,NA), momid=c(NA,"Mother",NA), 
                         sex = c("female", "female", "male"))
freqs = rep(1/6,6)
r = 0.003
model = 'equal'
als = c("15", "16", "16.2", "17", "18", "19")
locus1 <- FamiliasLocus(frequencies = freqs, 
                        name="Marker1", allelenames = als, 
                        MutationRate = r, MutationModel = model)
datamatrix <- data.frame(locus1.1=c("15","16.2","17"), locus1.2=c("18","18","19"))
rownames(datamatrix) <- persons
result <- FamiliasPosterior(ped2, locus1, datamatrix)
library(forrel)
#> Loading required package: pedtools
Familias2ped(ped2, datamatrix, locus1)
#> Error in pedmut::mutationModel(female = femalemut, male = malemut): unused arguments (female = femalemut, male = malemut)
Familias2ped(list(ped2), datamatrix, list(locus1))
#> Error in pedmut::mutationModel(female = femalemut, male = malemut): unused arguments (female = femalemut, male = malemut)

Created on 2018-09-26 by the reprex package (v0.2.0).

`readFam` and twins

In Windows Familas one can check if profiles from
tvilling1 and tvilling2 are identical by entering a parent child relationship and
overule by ticking 'Direct/identity' (!).
If the genotypes are identical, the likelihood is for one of them.
With a mutation model, the likelihood returned is 0.
readFam runs into problems:

library(forrel)
#> Loading required package: pedtools
readFam("http://familias.name/mz.fam")
#> Familias version: 3.2.8
#> Read DVI: No
#> 
#> Number of individuals (excluding 'extras'): 2
#>   Individual 'tvilling1': Genotypes for 2 markers read
#>   Individual 'tvilling2': Genotypes for 2 markers read
#> 
#> Number of pedigrees: 1
#> Warning in readFam("http://familias.name/mz.fam"): NAs introduced by coercion
#> Error in if (sex[par.idx] == 1) fidx.i[child.idx] = par.idx else midx.i[child.idx] = par.idx: missing value where TRUE/FALSE needed

^{Created on 2020-02-21 by the reprex package (v0.3.0)}

Missing error in `profileSim`

I suggest an error message is issued also in the last of the two below profileSim calls

library(forrel)
#> Loading required package: pedtools
x = nuclearPed(1)
profileSim(x, markers = marker(x, "1" = 1:2), id = 4)
#> Error: Unknown ID label: 4
profileSim(x, markers = marker(x), id = 4)
#> [[1]]
#>  id fid mid sex <1>
#>   1   *   *   1 -/-
#>   2   *   *   2 -/-
#>   3   1   2   1 -/-

^{Created on 2019-12-01 by the reprex package (v0.3.0)}

Error in missingPersonLR or profileSim in forrel v1.5.0?

Dear Magnus,

I was working with mispitools package, simLRgen function, and detected a bug when calculating LRs distributions considering an unrelated (poi1, lr1 in the function) and a related (poi2, lr2) poi. As you will be able to see, poi2 always gives low LR values, not consistent with what is expected (for example grandfather with 20 markers).

When I install a previous forrel version v1.4.1 everything works.
I have alerted about this issue in mispitools documentation.

Below I paste simLRgen function and a tiny example where you can run and detect the inconsistent result.

Bests,
Franco

simLRgen

simLRgen = function(reference, missing, numsims, seed) {
  st = base::Sys.time()

  if(pedtools::is.pedList(reference) && base::length(reference) == 1)
    reference = reference[[1]]

  if(!pedtools::is.ped(reference))
    base::stop("Expecting a connected pedigree as H1")

set.seed(seed)

poi1 = pedtools::singleton("poi1")
poi1 = pedtools::transferMarkers(from = reference, to = poi1)
poi1 = forrel::profileSim(poi1, numsims)

lr1 <- as.list(rep(NA, numsims))
for(i in 1:numsims) {
       lr1[[i]] = forrel::missingPersonLR(reference, missing, poi = poi1[[i]])
    }

poi2ped = forrel::profileSim(reference, numsims, ids = missing)

poi2 <- base::as.list(base::rep(NA, numsims))
for(i in 1:numsims) {
  poi2[[i]] = base::subset(poi2ped[[i]], missing)
}

base::rm(poi2ped)

lr2 <- base::as.list(rep(NA, numsims))

for(i in 1:numsims) {
  lr2[[i]] = forrel::missingPersonLR(reference, missing, poi = poi2[[i]])
}

LRsimulated <- base::cbind(base::sapply(lr1, function(x) {x[["LRtotal"]][["H1:H2"]]}),
   base::sapply(lr2, function(x) {x[["LRtotal"]][["H1:H2"]]}))
   base::colnames(LRsimulated) <- c("Unrelated", "Related")
   base::structure(base::as.data.frame(LRsimulated))
}

example

library(forrel)
x = linearPed(2)
plot(x)
x = setMarkers(x, locusAttributes = NorwegianFrequencies[1:5])
x = profileSim(x, N = 1, ids = 2)
datasim = simLRgen(x, missing = 5, 100, 123)

Unable to use know genotype data in pedigrees

I'm trying to get exclusionPower to run with known genotypes passed within the pedigrees. The following example fails, I think because the unrelated pedigree does not have marker data associated with it.

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))

m = marker(claim,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = c(0.5, 0.5))

claim = setMarkers(claim, m)

exclusionPower(claim, true, ids = c(1, 3), markerindex = 1, plot = FALSE)
#> Error: `partialmarker` must have length 1

^{Created on 2019-07-04 by the reprex package (v0.3.0)}

On the contrary, if I try the same example but directly supplying the genotypes as arguments to the exclusionPower function, I get the correct result:

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))

exclusionPower(claim, true,
               ids = c(1, 3),
               known_genotypes = list(c(3, 1, 1)),
               alleles = 2,
               afreq = c(0.5, 0.5),
               plot = FALSE)
#> [1] 0.25

^{Created on 2019-07-04 by the reprex package (v0.3.0)}

Table file input wizard

Consider creating a file wizard that lets the user choose what input format they want to provide.

Choose between tabs, commas or other characters as separators.
Integrate the module with table file inputs.
~~Allow multiple file input and concatenation.~~ Postponed until what other loading formats is clearly defined.
~~Choose which columns and/or rows do what in a file.~~ Allow users to transpose the table.
Predefine defaults for some of the most common input formats (e.g. for the relMix format for allele frecuencies)?
Possible heuristic for automatically choosing fields separators: the field separator is likely the most occurring whitespace character.

markerSim: Unknown ID label: mutmat

x = pedtools::nuclearPed(1)
forrel::markerSim(x, alleles = 1:2)
#> Unconditional simulation of 1 autosomal marker.
#> Individuals: 1, 2, 3
#> Allele frequencies:
#>    1   2
#>  0.5 0.5
#> Mutation model: No
#> Error: Unknown ID label: mutmat

Created on 2018-09-10 by the reprex package (v0.2.0).

readFam

In foo.fam read below, there is one marker and one ('unidentified') preson in the DVI module.
readFam return appears essentially empty:

library(forrel)
#> Loading required package: pedtools
x = readFam("http://familias.name/BookKEP/foo.fam",verbose = TRUE, useDVI =T)
#> Familias version: 3.2.8 (DVI)
#> Read DVI: Yes
#> 
#> Number of individuals (excluding 'extras'): 0
#> 
#> Number of pedigrees: 0
#> 
#> Database: New database
#> Number of loci: 1
#>   D3S1358: 9 alleles, model = equal, rate = 0 (unisex)
#> 
#> *** Reading DVI section ***
#> Returning the following families:
x
#> named list()

^{Created on 2020-02-09 by the reprex package (v0.3.0)}

exclusionPower could hang

This is probably obvious, but the documentation could
note that exclusionPower will hang if n is increased as indicated below:

library(forrel, quietly = TRUE)
n = 3 # hangs for large n, say n = 40
claim = nuclearPed(3, children = c("B1", "B2", "POI"))
x1 = nuclearPed(3, children = c("B1", "B2", "MP"))
x2 = singleton("POI")
true = list(x1,x2)
als = 1:n
p = rep(1/n, n)
exclusionPower(claim, true, ids = c("B1", "B2", "POI"), alleles = als, 
               afreq = p, plot = FALSE)
#> [1] 0.05144033

^{Created on 2018-10-22 by the reprex package (v0.2.1)}

Output from `transferMarkers`

transferMarkers should perhaps return an error below rather than returning
the last simulation?!:

library(forrel,quietly = TRUE)
v1 = singleton("V1")
m1 = marker(v1, alleles = 1:100, name = "L1")
v1 = setMarkers(v1, m1)
sim = profileSim(v1, N=2, cond = "L1", seed = 123)
transferMarkers(sim, v1, id = "V1")
#>  id fid mid sex    L1
#>  V1   *   *   1 42/90

^{Created on 2018-11-07 by the reprex package (v0.2.1)}

inbreeding and simulation

Something simillar (magnusdv/pedprobr#1) for the wish-list: simulation with inbred founders

library(pedtools)
library(forrel)
x = singleton(1)
m = marker(x, alleles= 1:5)
founder_inbreeding(x,1) = 1
markerSim(x, N=1 , partialmarker = m, seed = 17)
#> Unconditional simulation of 1 autosomal marker.
#> Individuals: 1
#> Allele frequencies:
#>    1   2   3   4   5
#>  0.2 0.2 0.2 0.2 0.2
#> Mutation model: No 
#> 
#> Simulation finished.
#> Number of calls to the likelihood function: 0.
#> Total time used: 0.05 seconds.
#>  id fid mid sex <1>
#>   1   *   *   1 2/1

Created on 2018-08-15 by the reprex package (v0.2.0).

kinshipLR

Very nice new version of kinshipLR()! One thing: apparently hypotheses can be given identical names and also be replicated, should that be allowed?

library(forrel)
#> Loading required package: pedtools
ids = c("A", "B")
sibs = nuclearPed(children = ids)
sibs = simpleSim(sibs, N = 1, ids = ids, seed = 123, alleles = 1:2, verbose =F)
unrel = list(singleton("A"), singleton("B"))
kinshipLR(H = sibs,  H = unrel)
#> Total LR:
#>  H:H  H:H 
#> 1.25 1.00
kinshipLR(H = sibs,  H = unrel, H = sibs, source = "H")
#> Total LR:
#> H:H H:H H:H 
#> 1.0 0.8 1.0

^{Created on 2020-05-20 by the reprex package (v0.3.0)}

markerSim

In the markerSim example
library(pedtools)
partial = marker(x, 3, 1, alleles=1:3)
Apparently the marker is created without genotype
Should it be ?:
partial = marker(x, "3"=1, alleles=1:3)

The below doesn't work unless mutmat is removed
library(forrel)
z = fullSibMating(2)
l = leaves(z)
z = addChildren(z, father = l[1], mother = l[2], sex = 2, nch=1)
p = c(1,1,1)/3; lp = length(p); alleler = 1:lp
R = 0.05
M = matrix(ncol= lp,nrow =lp, R/(lp-1), dimnames = list(alleler, alleler))
diag(M) = 1-R
m = marker(z, alleles = alleler, afreq = p, mutmat = M)
markerSim(z, N=1, partialmarker = m )

Custom pedigree not automatically selected when file is loaded

When selecting a file through the input field in the pedigrees tag, the user still has to choose the Custom option from the dropdown menu above.

Solution: make it so that when a pedigree file is loaded, the Custom option is selected automatically.

`Familias2ped`

I had some problems:

# Example 1. Converting Familias singletons to ped:
library(Familias, quietly = TRUE)
library(forrel, quietly = TRUE)
data(NorwegianFrequencies)
TH01 = NorwegianFrequencies$TH01
locus1 = FamiliasLocus(TH01)
persons = c('singleFemale')
datamatrix = data.frame(locus1.1 = 8, locus1.2= 9.3)
rownames(datamatrix) = persons
# First try
Familias2ped(NULL, datamatrix, locus1)
#> Error in data.frame(id = id, fid = 0, mid = 0, sex = 1, stringsAsFactors = FALSE): arguments imply differing number of rows: 0, 1
# Second try
ped1 = FamiliasPedigree(id = persons, dadid = NA,  momid = NA,  sex = 'female')
Familias2ped(ped1, datamatrix, locus1)
#> Error in matrix("0", nrow = nFath + nMoth, ncol = ncol(allelematrix)): non-numeric matrix extent

# Example 2. A small extension of the first
# example in the help for Familias2ped. Genotypes are lost when a singleton `NN` is added :

persons = c('mother', 'daughter', 'AF', 'NN')
ped1 = FamiliasPedigree(id = persons,
                        dadid = c(NA, 'AF', NA, NA),
                        momid = c(NA, 'mother', NA, NA),
                        sex = c('female', 'female', 'male', 'male'))
datamatrix = data.frame(THO1.1=c(NA, 8, NA, 8), THO1.2=c(NA,9.3, NA, 8))
rownames(datamatrix) = persons
x = Familias2ped(ped1, datamatrix, locus1)
plotPedList(x, marker = 1)

^{Created on 2019-01-11 by the reprex package (v0.2.1)}

Internal order and `markerSim`

Here's an error message I first ran into using profileSim

library(forrel, quietly = TRUE)
x = data.frame(id = 1:3, fid = c(3,0,0), mid = c(2,0,0),  sex = c(1, 2, 1))
x = as.ped(x)
m = marker(x, "3" = 1:2)
x = addMarkers(x, m)
markerSim(x, partialmarker = 1)
#> Note: Changing the internal order so that all parents precede their children.
#> Error: Input is not a `marker` object

This works with
x = data.frame(id = 1:3, fid = c(0,0,1), mid = c(0,0,2), sex = c(1, 2, 1))
^{Created on 2019-01-12 by the reprex package (v0.2.1)}

ibd_identity problem

library(ribd)
library(pedtools)
x = halfSibStack(2)
ibd_identity(x, ids = c(4,5))
#> Error in `[<-`(`*tmp*`, i, i, i, value = (1 + 3 * fi)/4): subscript out of bounds

Created on 2018-08-13 by the reprex package (v0.2.0).

kinshipLR problem

I have a problem with kinshipLR. Here's the first example of the documentation:

library(forrel)
#> Loading required package: pedtools
ids = c("A", "B")
sibs = nuclearPed(children = ids)
sibs = simpleSim(sibs, N = 5, alleles = 1:4, ids = ids, 
                 seed = 123, verbose = FALSE)
halfsibs = relabel(halfSibPed(), old = 4:5, new = ids)
unrel = list(singleton("A"), singleton("B"))
kinshipLR(sibs, halfsibs, unrel)
#> Error in vapply(markers, function(i) likelihood(xx, marker1 = i), FUN.VALUE = 1): values must be length 1,
#>  but FUN(X[[1]]) result is length 5

^{Created on 2020-07-06 by the reprex package (v0.3.0)}

Familias2ped - `please contact MDV

Below, based on first example in the documentation of
Familias2ped, I am adviced to contact MDV and I hereby comply:

library(Familias, quietly = TRUE)
library(forrel)
#> Loading required package: pedtools
data(NorwegianFrequencies)
TH01 = NorwegianFrequencies$TH01
locus1 = FamiliasLocus(TH01)
persons = c('mother', 'daughter', 'AF')
ped1 = FamiliasPedigree(id = persons,
                        dadid = c(NA, 'AF', NA),
                        momid = c(NA, 'mother', NA),
                        sex = c('female', 'female', 'male'))
datamatrix = data.frame(THO1.1=c(NA, 8, NA), THO1.2=c(NA,9.3, NA))
rownames(datamatrix) = persons
x = Familias2ped(ped1, datamatrix, locus1)
#> Warning in as.ped.data.frame(p, locus_annotations = annotations): Argument
#> `locus_annotations` is deprecated; use `locusAttributes` instead
#> Error: Genotype columns are sorted differently from `locusAttributes`. Please contact MDV
```

<sup>Created on 2019-08-11 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)</sup>

readFam problem

After quite a bit of testing I manage to to find a problem:

x = forrel::readFam("http://familias.name/foo.fam")
#> Familias version: 3.2.9 (DVI)
#> Read DVI: Yes
#> 
#> Number of individuals (excluding 'extras'): 0
#> 
#> Number of pedigrees: 0
#> 
#> Database: New database
#> Number of loci: 1
#>   D3S1358: 9 alleles, model = equal, rate = 0 (unisex)
#> 
#> *** Reading DVI section ***
#> Unidentified persons: 2
#>   BS1
#>   BS2
#> 
#> Reference families: 1
#>   F14 (3 persons, 2 pedigrees)
#>     Reference pedigree
#> Error in this.fidx[child.idx[par.is.male]] <- parent.idx[par.is.male]: NAs are not allowed in subscripted assignments

^{Created on 2020-02-12 by the reprex package (v0.3.0)}

kinshipLR complains about missing linkage map

kinshipLR complains about missing linkage map for a case with one marker.
The problem goes away if the unnecessary call to transferMarkers is dropped.
library(forrel)
#> Loading required package: pedtools
library(ibdsim2)
x = nuclearPed(1)
pattern = ibdsim(x, 1, verbose = F)
x = setSNPs(x, forrel::FORCE[1:2,])
id = c(1,3)
H1 = profileSimIBD(x, pattern, ids = id, seed = 123)[[1]]
H2 = list(singleton(id[1]), singleton(id[2]))
H2 = transferMarkers(H1, H2, ids = id)
kinshipLR(H1, H2, markers = 1)
#> Error: Linked markers detected, but no linkageMap provided
Created on 2022-03-29 by the reprex package (v2.0.1)

LRpower output

I believe that there is a typo in the last line of the output below and that it should be
P(LR ≥ 0)

library(forrel)
#> Loading required package: pedtools
PO = nuclearPed(fa = "A", ch = "B")
UN = list(singleton("A"), singleton("B"))
LRpower(PO, UN, UN, ids = c("A", "B"), nsim = 1, threshold = 0, 
        alleles = 1:100, seed = 123, verbose = F)
#> LR power analysis using individuals A, B.
#> Mean total LR: 0 
#> Mean total log10(LR): -Inf 
#> Estimated inclusion power: 
#>   P(LR > 0) = 1

^{Created on 2020-08-26 by the reprex package (v0.3.0)}

Marriage loops and `markerSim`

The help file could note that marriage loops are not implemented in full
generality in markerSim and hence profileSim.

library(forrel)
#> Loading required package: pedtools
library(pedprobr)
x = nuclearPed(2, sex = c(2,2))
x= addChildren(x, mother = 3, nch = 2)
#> Father: Creating new individual with ID = 5
x= addChildren(x, mother = 4,nch = 1,father = 5)
m = marker(x, "4" = 1:2)
x = setMarkers(x, m)
markerSim(x, ids = 3, partialmarker = 1)
#> Conditional simulation of 1 autosomal marker.
#> Target individuals: 3
#> Conditioning on the following data:
#>    <NA>
#>  1  -/-
#>  2  -/-
#>  3  -/-
#>  4  1/2
#>  5  -/-
#>  6  -/-
#>  7  -/-
#>  8  -/-
#> -------- 
#> Chrom NA: NA (Mb), NA (cM)
#> Mutation model: No 
#> Allele frequencies:
#>    1   2
#>  0.5 0.5
#> Marriage loops detected, trying different selection method
#> Loop breakers: 4
#> No further loops to break
#> 
#> Simulation strategy
#> ===================
#> Pre-computed joint distribution: 3
#> Brute force conditional simulation: None
#> Hardy-Weinberg sampling (founders): None
#> Simple gene dropping: None
#> Required likelihood computations: 3
#> Error in attributes(mk) <- attrib: dims [product 18] do not match the length of object [16]
profileSim(x, ids = 3, partialmarker = 1)
#> Error in markerSim(x, N = N, ids = ids, partialmarker = pm, verbose = F, : formal argument "partialmarker" matched by multiple actual arguments

x = breakLoops(x, 3)
#> Loop breakers: 3
markerSim(x, ids = 3, partialmarker = 1)
#> Error: `ped` objects with pre-broken loops are not allowed as input to `markerSim()`
profileSim(x, ids = 3, partialmarker = 1)
#> Error in markerSim(x, N = N, ids = ids, partialmarker = pm, verbose = F, : formal argument "partialmarker" matched by multiple actual arguments

^{Created on 2019-12-26 by the reprex package (v0.3.0)}

Simulation output

Below is an example where likelihood calculations are compared to simulations.
For me the comparison would have been simpler if genotypes from
simpleSim and markerSim had been sorted so, i.e not giving both 1/2 and 2/1.
Or is there is a simple way to deal this? gsub below seems clumsy:

library(pedtools)
library(pedprobr)
library(forrel)
p = c(0.5, 0.5); lp = length(p); alleler = 1:lp
Rfe = 0.05; Rma = 0.05 
Mfe = mutationMatrix(alleler, "eq", Rfe)
Mma = mutationMatrix(alleler, "eq", Rma)
z = fullSibMating(1)
l = leaves(z)
z = addChildren(z, father = l[1], mother = l[2], sex = 2, nch=1)
m = marker(z, alleles = alleler, afreq = p, mutmat = list(female = Mfe, male = Mma))
chrom(m) = "X"
l = leaves(z)
one = oneMarkerDistribution(z, m, ids = l)
Nsim = 10000
sim1 = markerSim(z, N = Nsim, partialmarker = m, available = l )
# sim1 = simpleSim(z, Nsim, alleler, p,5)
foo = as.data.frame(sim1)[l, -(1:4)]
foo = gsub("2/1", "1/2", foo)
one.sim = table(foo)/Nsim
one[c(1,3,2)] - one.sim

Single marker frequency database files do not work

This is a strange error. When loading a tabular file as the frequency database, if that file happens to contain only frequency information for one marker, exclusionPower calculation fails.

This is an example of an offending file:

"","M1"
"1",0.5
"2",0.5

More specifically, what files is obtaining the allele denominations of that marker, which is done by computing

alleles = rownames(frequencyDB()[!is.na(frequencyDB()[markerName]),])

Where markerName is "M1" and frequencyDB() is a well-formed dataframe.

A file with two markers does not cause this warning.

IBDtriangle

Change IBDtriangle to allow for k0 and k2 on the axes rather than the current greek kappas, eg by changing the argument to allow: kappas = c("k0", "k2)?

Keep marker names in `profileSim`

Should profileSim keep marker names as opposed to what's currently implemented?:

library(forrel, quietly =TRUE)
x = nuclearPed(children = c("B1", "B2"))
m1 = marker(x, B1 = 1:2, alleles = 1:3, afreq = c(.2, .3, .5), name ="locus1")
profileSim(x, N = 1, ids = "B2", conditions = list(m1))
#> [[1]]
#>  id fid mid sex <1>
#>   1   *   *   1 -/-
#>   2   *   *   2 -/-
#>  B1   1   2   1 1/2
#>  B2   1   2   1 2/3

^{Created on 2018-11-05 by the reprex package (v0.2.1)}

Rearranging `profileSim` output

Below I explore the profileSim function and rearrange the output allowing for transferMarkers
and likelihood calculations. Any comments on improved code is appreciated.

library(forrel,quietly = TRUE)
fa = singleton("F")
m = marker(fa, name = "L1")
fa = setMarkers(fa, m)
ch = singleton("C")
m = marker(ch,  name = "L1")
ch = setMarkers(ch, m)
nsim = 5
set.seed(123)
sim = profileSim(list(fa, ch), N = nsim, ids = c("F", "C"), cond = "L1")
simN = vector("list", nsim)
for (i in 1:nsim)
  simN[[i]] = lapply(sim, function(z) z[[i]])
LR(simN, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3]   [,4] [,5]
#> L1 0.25 0.125 0.25 0.0625 0.25
x = nuclearPed(1, father = "F", children = "C")
m = marker(x,  name = "L1")
x = setMarkers(x,m)
for (i in 1:nsim)
  simN[[i]] = transferMarkers(simN[[i]], x, id = c("F","C"))
LR(simN, 1)$likelihoodsPerSystem
#>    [,1]  [,2] [,3] [,4] [,5]
#> L1 0.25 0.125 0.25    0 0.25

^{Created on 2018-11-07 by the reprex package (v0.2.1)}

EP: mutation only possible toobserved alleles

The code below returns EP = 0; I thought EP > 0

library(forrel)
#> Loading required package: pedtools
z = nuclearPed(4, children = c(3, 4, 5, "POI"))
p = c(0.1, 0.2, 0.3, 0.4)
als = 1:4
m = marker(z, alleles = als, afreq = p, "3" = 1, "4" = 2, "5" = c(1,3) )
z = addMarkers(z,m)
#plot(z, 1, skip.empty.genotypes = TRUE)
r = 0.06
M = rbind(c(1-r, r/2, r/2, 0),
          c(r/2, 1-r, r/2, 0),
          c(r/2, r/2, 1-r, 0),
          c(0, 0,  0, 1))
dimnames(M) = list(1:4, 1:4)
mutmod(z,1) = list("custom", matrix = M)
missingPersonEP(z, "POI", disableMutations = FALSE, verbose = FALSE)
#> $EPperMarker
#> 1 
#> 0 
#> 
#> $EPtotal
#> [1] 0
#> 
#> $expectedMismatch
#> [1] 0
#> 
#> $distribMismatch
#> 0 
#> 1
# Thore thinks EP should be:
p[4]^2+2*p[4]*(1-p[4])
#> [1] 0.64

^{Created on 2019-10-30 by the reprex package (v0.3.0)}

'ids' argument in showInTriangle

Assume that IBD coefficients have been estimated between all pedigree members. The output from e.g. IBDestimate() can then be used directly as input to showInTriangle() and the set of coefficients for each pair are plotted.

My question is then: Is it possible to specify a subset of pedigree member pairs to plot when using showInTriangle()? I cannot find 'ids' as an option in showInTriangle. This feature would be helpful.

x = nuclearPed(father = "fa",mother = "mo", children = c("ch1","ch2"), sex = c(1,2))
mlist = lapply(1:4, function(i) marker(x, alleles = letters[1:10], name = paste0("M",i)))
x = setMarkers(x, mlist)
x = profileSim(x,N = 1)
k = IBDestimate(x)
IBDtriangle()

Suggestion of code:

Only want to plot coefficients between ch1 and ch2:
showInTriangle(k, ids = c("ch1","ch2"))

Or if more pairs are wanted, a matrix of ids may be specified:
showInTriangle(k, ids = rbind(c("ch1","ch2"),c("fa","mo")))

`markerindex` problem in `exclusionPower

I have problems selecting the marker in exclusionPower;
it seems that markerindex = 1 always, also in plot (great, not shown in reprex) rendered by
exclusionPower:

library(forrel, quietly =T)
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))
p1 = c(0.5, 0.5)
m1 = marker(claim,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = p1)
claim = addMarkers(claim, m1)
p2 = c(0.25, 0.75)
m2 = marker(claim,
            `3` = c(2, 2),
            alleles = c(1, 2),
            afreq = p2)
claim = addMarkers(claim, m2)
PE1 = exclusionPower(claim, true, ids = 1, markerindex = 1,
                     verbose = FALSE)
PE2 = exclusionPower(claim, true, ids = 1, markerindex = 2,
                     verbose = FALSE)
PE1 == PE2 # should be false
#> [1] TRUE
PE2 == p2[1]^2 # should be TRUE
#> [1] FALSE

^{Created on 2019-07-07 by the reprex package (v0.3.0)}

Read frequency database on Familias and FamLinkX format

It would be nice to be able read tab separated asci files exported from Windows Familias and FamLinkX on the below format to define ped markers

DXS10148
13.3 0.4
14 0.3
17 0.3

DXS10135
16 0.5
16.1 0.5
...

Warnings from shownInTriangle

The plot below is fine, but I get some warnings,
maybe something local to my computer.

library(forrel)
#> Loading required package: pedtools
showInTriangle(k0 = 3/8, k2 = 1/8, label = "3/4 siblings", pos = 1)
#> Warning in par(opar): graphical parameter "cin" cannot be set
#> Warning in par(opar): graphical parameter "cra" cannot be set
#> Warning in par(opar): graphical parameter "csi" cannot be set
#> Warning in par(opar): graphical parameter "cxy" cannot be set
#> Warning in par(opar): graphical parameter "din" cannot be set
#> Warning in par(opar): graphical parameter "page" cannot be set

^{Created on 2020-01-14 by the reprex package (v0.3.0)}

Error when plotting pedigrees with `exclusionPower()`

Hi!

I ran into this problem when trying to run exclusionPower. Below is the output from reprex.

library(pedtools)
library(paramlink)
#> 
#> Attaching package: 'paramlink'
#> The following objects are masked from 'package:pedtools':
#> 
#>     addDaughter, addParents, addSon, ancestors, branch,
#>     breakLoops, connectedComponents, cousinPed, cousins,
#>     descendants, doubleCousins, doubleFirstCousins,
#>     findLoopBreakers, findLoopBreakers2, fullSibMating,
#>     getMarkers, grandparents, halfCousinPed, halfSibStack,
#>     is.singleton, leaves, marker, mendelianCheck, mergePed,
#>     nephews_nieces, nuclearPed, offspring, parents, plotPedList,
#>     quadHalfFirstCousins, randomPed, relabel, removeIndividuals,
#>     removeMarkers, setMarkers, siblings, singleton, spouses,
#>     swapSex, tieLoops, unrelated
library(forrel, quietly = T)
#> 
#> Attaching package: 'forrel'
#> The following objects are masked from 'package:paramlink':
#> 
#>     IBDestimate, IBDtriangle, LR, exclusionPower, markerSim,
#>     readFamiliasLoci, showInTriangle, simpleSim

############################################
### A standard case paternity case:
### Compute the power of exclusion when the claimed father is in fact
### unrelated to the child.
############################################

# Claim: Individual 1 is the father of indiv 3
claim = nuclearPed(nch = 1, sex = 2)
#> Error in nuclearPed(nch = 1, sex = 2): unused argument (nch = 1)

# Truth: 1 and 3 are unrelated
true = list(singleton(id = 1), singleton(id = 3, sex = 2))

# Indivs 1 and 3 are available for genotyping
ids = c(1, 3)

# An equifrequent SNP
als = 2
afreq = c(0.5, 0.5)

# The exclusion power assuming no known genotypes
PE1 = exclusionPower(claim, true, ids = ids, alleles = als, afreq = afreq)
#> Error in is.ped(ped_claim): object 'claim' not found

^{Created on 2019-06-05 by the reprex package (v0.3.0)}

However, if I don't load the pedtools and paramlink packages I get a different error:

library(forrel, quietly = T)

############################################
### A standard case paternity case:
### Compute the power of exclusion when the claimed father is in fact
### unrelated to the child.
############################################

# Claim: Individual 1 is the father of indiv 3
claim = nuclearPed(nch = 1, sex = 2)

# Truth: 1 and 3 are unrelated
true = list(singleton(id = 1), singleton(id = 3, sex = 2))

# Indivs 1 and 3 are available for genotyping
ids = c(1, 3)

# An equifrequent SNP
als = 2
afreq = c(0.5, 0.5)

# The exclusion power assuming no known genotypes
PE1 = exclusionPower(claim, true, ids = ids, alleles = als, afreq = afreq)
#> Error in kinship2::pedigree(id = ped$id, dadid = ped$fid, momid = ped$mid, : Wrong length for affected

^{Created on 2019-06-05 by the reprex package (v0.3.0)}

I was using R version 3.5.0 (2018-04-23) -- "Joy in Playing". As for the packages, I was using these versions:

> packageVersion('paramlink')
[1] '1.1.2'
> packageVersion('pedtools')
[1] '0.1.0'
> packageVersion('forrel')
[1] '0.1.0'

Thanks for your time!

`markerSim()` error in inbred pedigree

library(forrel, quietly = T)

x = cousinPed(0, child = T)
x = setMarkers(x, marker(x, "3" = 1:2))

markerSim(x, partial = 1, verbose = F)
#> Error: Could not determine sensible order for gene dropping.

return from `getAlleles`

# `getAlleles` reportedly returns a
# `character matrix`, this is not the case below.
library(forrel, quietly = T)
v1 = singleton("V1")
m = marker(v1, "V1" = 1:2)
v1 = setMarkers(v1,m)
v2 = singleton("V2")
m = marker(v2, "V2" = 1:2)
v2 = setMarkers(v2,m)
from = list(v1, v2)
g = getAlleles(from, ids = "V1")
is.matrix(g)
#> [1] FALSE

^{Created on 2019-02-27 by the reprex package (v0.2.1)}

Description of profileSim

Something is wrong in the Description of this function:
" .... which does the is function is a"

Genotypes are lost in conversion from Familias

Moved from #13 (and slightly edited):

A small extension of the first example in the help for Familias2ped. Genotypes are lost when a singleton NN is added :

library(forrel, quietly = TRUE)
library(Familias, quietly = TRUE)

persons = c('mother', 'daughter', 'AF', 'NN')

ped1 = FamiliasPedigree(id = persons,
                        dadid = c(NA, 'AF', NA, NA),
                        momid = c(NA, 'mother', NA, NA),
                        sex = c('female', 'female', 'male', 'male'))

datamatrix = data.frame(THO1.1 = c(NA, 8,   NA, 8), 
                        THO1.2 = c(NA, 9.3, NA, 8),
                        row.names = persons)

# Genotypes are lost in conversion:
Familias2ped(ped1, datamatrix, loci = NULL)
#> $`_comp1`
#>        id fid    mid sex <1>
#>    mother   *      *   2 -/-
#>  daughter  AF mother   2 -/-
#>        AF   *      *   1 -/-
#> 
#> $`_comp2`
#>  id fid mid sex <1>
#>  NN   *   *   1 -/-

^{Created on 2019-01-16 by the reprex package (v0.2.1)}

magnusdv / forrel Goto Github PK

forrel's People

Contributors

Stargazers

Watchers

forrel's Issues

References

simLRgen

example

Recommend Projects

Recommend Topics

Recommend Org