Introduction to PALM

PALM is a software that conducts gene-based association analysis. PALM exploits Hidden Markov Models (HMM) to capture the inner gnomic structures. Then, unique gene features (termed "palm prints") are extracted for each gene set by calculating the natural gradient. Then palm-prints are tested with $\chi^2$ significance test.

PALM toolset is flexible in the gene data format. The input can be a entire data matrix, or spliting big data to small parts for the sake of limited RAM. Only small changes in the code are needed. Here, we assume a genome scale data need to be processed. Therefore, it can be hardly to load all the data at once. We split big dataset to small files and only load one file at each time.

Herein, we prepared a small dataset extract from CEU data for demonstration of PALM. This dataset includes 579 genes.

The main procedure for gene-based analysis includes:

• Data preparation
• Gene set assembling
• HMM training
• Palm-prints extracting
• Significance testing

To use the software, we first need to include the function files in the MATLAB path. Use:

addpath('./PALM_toolset/','./PALM_toolset/HMM_functions/');

We will demonstrate the details step by step. All the commands are in the demo.m file.

Date preparation

In this demonstration, we conduct gene-based association analysis based on a small dataset extracted from Utah Residents with Northern and Western European ancestry cohort (CEU) in the 1000 Genomes Project (CEU_sub). This dataset spans 579 gene regions.

For phenotype, we follow the instructions in our paper and simulate case-control statuses based on a regression model.

$logit(y=1) = \alpha x + \beta c_i + \varepsilon$

Here, we randomly select one SNP in the $20^{th}$ gene set as causal marker and generate simulated phenotype.

Genotype file

PLINK format is a widely used form for storing gene data (either in .ped or .bed). Here we assume the original data were stored in this format. Firstly, we recode the data in additive model using PLINK by command:

plink --file CEU_sub --recodeA --out CEU_sub

Here we assume the PLINK has been in the system path. For more information, please refer to the PLINK website. Then, a file named CEU_sub.raw is obtained. The successive process is based on the additive genotype.

We have implemented this command into a Windows batch file (recode_raw.bat) for users' convenience. Users only need to replace the file name (CEU_sub) to their owns.

Phenotype file

Phenotype is generally deposited in .bim file or .bed file. Here we simply store it in a plain text file (pheno.txt). The format is a $ n\times 1$ vector where $n$ is sample number. This can be easily achieved using PLINK.

Covariants

PALM conducts significance test using a regression based model. Therefore, it can be convenient to incorporate covariants for analysis. The covariants data also deposit in matrix using plain text format with shape $n \times m$ where $n$ again is a sample size and $m$ is the number of covariants. In the demo, we only consider gender which is recorded in file gender.txt

Gene set assembling criterion

PALM can adapt to different criteria for grouping gene sets as long as following the sequence orders. We use the gene location as way to assemble sets. The RefGene (hg19) information is obtained from UCSC database. We have processed it into a MATLAB format (gene_location.mat) as well as a simple text format (gene_location.txt).

For users who want to use their own set partition criterion, they can make a set file using MATLAB .mat format. The gene set file is a $n\times 3$ MATLAB cells where $n$ is the number of sets. The first column is the set index (just for index and will not use in the the later analysis), and the second and third columns are set content and chromosome index. In the second column, there is a vector containing all the SNP indexes belonging to this set. The SNP index is the row number in the PLINK-formatted SNP information file .map.

Gene set assembling

We provide two file for gene set assembling: SNP2gene.m and assemble_gene_set.m. Generally, these two files read the gene information from gene location file and obtain SNP location from CEU_sub.map. Then, generates the corresponding gene sets based on the given location. To take the regulatory impacts into consideration, each data set also contained the down-stream and up-stream SNPs within 5 kb of a gene.

To encode the gene set file, use

assemble_gene_set(gene_location_path, map_path, gene_set_name, set_filter)

where

gene_location_path: The file path and name of the gene location file. map_path: The path of the standard pink *.map file of the test dataset. gene_set_name: Indicating the output path and file name of the assembled gene set file. set_filter: A flag to assign the minimum set size.

Please refer to the function file for details of the function parameters. SNP2gene.m works in a similar way except it works on a parameter-free mode. Details can also be found in its function file.

In the demo.m, the commands as following are conducting this process:

gene_location_path='gene_location.txt'; map_path='CEU_sub.map'; gene_set_name='CEU_set.mat'; set_filter=1;

assemble_gene_set(gene_location_path,map_path,gene_set_name,set_filter)

HMM training

With gene sets, we are able to train an HMM based on each set. To speed the training process, we here adopt a split-train strategy that training several HMM models on different data in parallel with the same initiation. To do this, it would be better to deposit each gene set in an independent file. This function is completed using recode_data_split() function.

recode_data_split(raw_path,geno_path,gene_set_name,individual_num)

where

raw_path: The path of the CEU_sub.raw. geno_path: The destination of generated .geno files. gene_set_name: The path of gene set file (gene_set.mat). individual_num: Sample numbers in each file.

In the demo.m, this process is implemented by:

raw_path='CEU_sub.raw'; geno_path='./CEU_geno/CEU_geno'; individual_num=200;

recode_data_split(raw_path,geno_path,gene_set_name,individual_num)

Please refer to the comments in recode_data_split.m file for details.

After generated .geno files, we input the gene data to HMM for training:

O=3; Q=5; HMM_para_path='CEU_hmm.mat';

train_HMM(O,Q,geno_path,individual_num,HMM_para_path,gene_set_name,'Training interval',30,'Parallel','Enable')

In the function train_HMM():

O: Number of genotype category. Q: Number of hidden nodes in HMM. geno_path: The output path of .geno file by recode_data_split. individual_number: Number of samples in the training dataset. HMM_para_path: The aim path and .mat name of the HMM parameters. gene_set_name: The gene_set path generated from assemble_gene_set function.

Keyword 'Gene frequncy': if missing genotypes exist in the data. Keyword 'Training interval': if Training interval is given, we can split all the data to several groups and training them at the same time. Keyword 'Parallel': whether we need start multiple CPU kernels in parallel in MATLAB

Users can find details in the file train_HMM.m.

Palm-prints extracting

In the training process, we obtain the HMM parameters deposited in HMM_para_path. We next compute the natural gradient mapping and extract "palm-prints". Related codes are:

HMM_para.prior = prior; HMM_para.transmat = transmat; HMM_para.obsmat = obsmat;

extract_PALM_print(HMM_para,geno_path, gene_set, individual_num, PALM_path)

where

HMM_para: Parameters obtained form HMM training, including initial probability $Q\times 1$; emission matrix $Q\times O$; transition matrix $Q\times Q$ where $Q$ is the hidden states number, $O$ is the category of observed states. geno_path: The output path of .geno file by recode_data_split. gene_set: The gene_set file specify the how SNPs are grouped. individual_sum: Number of samples in the training dataset. PALM_path: The desired output path of palm-print. Each set is in an independent .mat file.

Users can refer to extract_PALM_print.m for details.

Significance testing

With the obtained palm-prints, we can further conduct association test with phenotype. The significance test is the conventional regression-based test. In demo.m, we calculate the $P$ value for each set in a loop manner:

for set_count=1:gene_number select_PALM_print=importdata([PALM_path,'_',... num2str(set_count),'.txt']); p_val=association_test(phenotype,X,select_PALM_print); result_p_value(set_count,:)=p_val; disp(['Computing set ',num2str(set_count)]); end

where

phenotype: A $n\times 1$ phenotype file coded in 0-1 where 1 indicates the disease. $n$ is the number of individuals. cov: A $n\times c$ matrix recording the covariants. $c$ is the number of factors need to be adjusted. palm_print: The palm_print need to be tested. $n\times p$ matrix where $p$ is the dimension of the palm-print.

The output $P$ values for all sets are located in result_p_value variable.

Application to other datasets

For readers who want to run their own datasets, they can directly modify the demo.m file by replacing the aim path in the file to their own. Section Date preparation has presented the detail information, authors need to modify:

Genotype files

Phenotype files

Covariants matrix

Gene set allocation