bbitarello / ncd-statistics Goto Github PK

Author: Bárbara D. Bitarelo

Created: 27.03.2015

Last modified: 28 June 2022

Language: R

Welcome to the NCD repo!

UPDATE (June 2022): We have an R package under development that implements NCD1 and NCD2. Check it out here:

This repository provided scripts related to the article "Signatures of long-term balancing selection in human genomes":

What are the NCD (Non-Central Deviation) statistics?

NCD statistics measure the average difference between allele frequencies in a given region from a deviation point, which we call the 'target frequency (tf)'. So, assuming tf = 0.5, the more the allele frequencies are close to 0.5, the lower the NCD values. We propose two implementations of this statistic: NCD1 (only SNPs are required and used as informative sites) and NCD2 (SNPs and fixed differences are used as informative sites). In the manuscript, NCD2 as used to scan the human genome.

Here, we show how to:

• item Run NCD1 and NCD2
• item Using examples from the manuscript (see above) which can be extended to other species

Running NCD:

this requires:

* SNP input data as a data.table object in R

* Fixed differences (FD) input data (e.g. human-chimp FD bed file, as used in the manuscript) as a data.table object in R.
Note: *NCD1* only requires the first input file, whereas NCD2 requires both
    
* An open R session

Example input files are provided in example_input_files/. Refer to the README.md in that directory for further explanations.

First:

• clone this repo: go to your directory and clone:

git clone https://github.com/bbitarello/NCD-statistics.git

Note: this may take a few minutes because large example files are provided!

• go to root NCD directory

cd NCD-statistics/

Second:

• open an R session and type

source('scripts/preamble.R') #loads several packages
source('scripts/NCD_func.R') #loads NCD functions NCD1 and NCD2
readRDS('example_input_files/SNP_test_input.rds')-> SNP_input #necessary for NCD1 and NCD2
readRDS('example_input_files/FD_test_input.rds')-> FD_input  #only necessary for NCD2
system.time(example.run.ncd1<-foreach(x=1:22, .combine="rbind", .packages=c("data.table")) %dopar% NCD1(X=SNP_input[[x]], W=3000, S=1500)); #  6 seconds 
system.time(example.run.ncd2<-foreach(x=1:22, .combine="rbind", .packages=c("data.table")) %dopar%  NCD2(X=SNP_input[[x]], Y=FD_input[[x]],  W=3000, S=1500)); # 9 seconds

Note that the runtime will vary considerably depending on computational constraints. My experience with registerDoMC(11) is that each scan for the entire genome takes about one minute. The examples here encompass a smaller proportion of the genome and run in a few seconds. See example_input_files/README.md

Important: This is an example. It is a roadmap of how NCD can be used for other input data (including non-SNP data). Even though the example input data is real (phase 3 1000 Genomes), it does not reproduce the findings from the paper.

Acnowledgements

Many thanks to @VitorAguiar (https://github.com/VitorAguiar) with optimization of the NCD codes.

Updates

Soon: will fix NCD2 so that it runs even when there are no FDs in the input file.(update: fixed) Update: Fixed issue with NCD2 output.