Test for association between a trait and principal components of genotypes within a region, on summary statistics
PCA(scoreFile, geneFile, regions, cor.path = "", annoType = "",
n, beta.par = c(1, 1), weights.function = ifelse(maf > 0,
dbeta(maf, beta.par[1], beta.par[2]), 0), var.fraction = .85,
write.file = FALSE)
name of data file generated by prep.score.files()
.
name of a text file listing genes in refFlat format. If not set, hg19 file will be used (see Examples below).
character vector of gene names to be analysed. If not set, function will
attempt to analyse all genes listed in geneFile
.
path to a folder with correlation files (one file per each gene to be analysed).
Names of correlation files should be constructed as "geneName.cor" (e.g. "ABCG1.cor", "ADAMTS1.cor", etc.)
Each file should contain a square matrix with correlation coefficients (r) between genetic variants
of a gene. An example of correlation file format:
"snpname1" "snpname2" "snpname3" ...
"snpname1" 1 0.018 -0.003 ...
"snpname2" 0.018 1 0.081 ...
"snpname3" -0.003 0.081 1 ...
...
One way to generate such file from original genotypes is:
write.table(cor(g), file = paste0(geneName, ".cor"))
where g
is a genotype matrix (nsample x nvariants) for a given gene with genotypes coded as 0, 1, 2
(exactly the same coding that was used to generate betas).
for files annotated with the seqminer
package, a character (or character vector) indicating annotation types to be used (e.g.
"Nonsynonymous", "Start_Loss", "Stop_loss", "Essential_Splice_Site")
size of the sample on which summary statitics were obtained.
two positive numeric shape parameters in the beta distribution to assign weights for each genetic variant as a function of MAF in the default weights function (see Details). Default = c(1, 1).
a function of minor allele frequency (MAF) to assign weights for each genetic variant. By default, the weights will be calculated using the beta distribution (see Details).
minimal proportion of genetic variance within region that should be explained by principal components used (see Details for more info).
output file name. If specified, output (as it proceeds) will be written to the file.
a data frame containing P values, numbers of variants and filtered variants for each of analyzed regions. It also contains the number of the principal components used for each region and the proportion of genetic variance they make up.
PCA test is a useful tool to detect association between genetic variants of a region and a trait
when genetic variants are strongly correlated. PCA test is based on the spectral decomposition of
correlation matrix among genetic variants. The number of top principal components will be chosen
in such a way that >= var.fraction
of region variance can be explained by these PCs.
By default, var.fraction
= 0.85, i.e. >= 85% of region variance will be explained by PCs.
If var.fraction
= 1 then the results of PCA test and MLR-based test are identical.
beta.par = c(a, b)
can be used to set weights for genetic variants.
Given the shape parameters of the beta function, beta.par = c(a, b)
,
the weights are defined using probability density function of the beta distribution:
\(W_{i}=(B(a,b))^{^{-1}}MAF_{i}^{a-1}(1-MAF_{i})^{b-1} \),
where \(MAF_{i}\) is a minor allelic frequency for the \(i^{th}\) genetic variant in the region, which is estimated from genotypes, and \(B(a,b)\) is the beta function. This way of defining weights is the same as in original SKAT (see [Wu, et al., 2011] for details). Precision of input values (betas and P) can be important for perfect correspondence between PCA on summary statistics and PCA performed on genotypes. We suggest to keep as much decimal points as possible for input values.
Jolliffe, I.T. A note on the use of principal components in regression. J R Stat Soc Ser C 31, 300-303 (1982).
# NOT RUN {
## Run PCA with example files:
VCFfileName <- system.file("testfiles/CFH.scores.anno.vcf.gz",
package = "sumFREGAT")
cor.path <- system.file("testfiles/", package = "sumFREGAT")
n <- 85 # your sample size
out <- PCA(VCFfileName, region = 'CFH', cor.path = cor.path, n = n)
# }
Run the code above in your browser using DataLab