Population-wise decay of linkage disequilibrium (LD)
0
0
Entering edit mode
8.1 years ago
akang ▴ 110

I am trying to calculate the extent of LD in 10 populations. I have a plink file containing all the data in one ped file. Below is the command i used to calculate the ld and then plot the results in R. But the plot shows everything as one. I dont get the population wise plot.I am certainly doing alot wrong here. Do i need to specify family id here? How to go about it?Ill appreciate any help.

--noweb
--ped Input.txt
--map Input.map
--no-fid
--no-parents
--no-sex
--no-pheno
--compound-genotypes
--missing-genotype 9
--allow-no-sex
--maf 0.01
--hwe 0.0001
--geno 0.1
--mind 0.1
--cow
--r2
--ld-window 99999
--ld-window-r2 0
R
ld <- read.table("plink.ld",header=T)
a<-(ld$BP_B-ld$BP_A)
plot(a,ld$R2)
Plink plot LD • 4.1k views
ADD COMMENT
0
Entering edit mode

May be you can solve it on the level of R... How does the plink.ld file looks like?

ADD REPLY
0
Entering edit mode

It has 7 columns like below and 175288 rows

    CHR_A   BP_A           SNP_A CHR_B   BP_B           SNP_B          R2
   1 244698 SNP1     1 524506 SNP2 0.026080900
   1 244698 SNP1     1 643558 SNP3 0.029830100
ADD REPLY
0
Entering edit mode

It is not too helpful. Anyway, if you have some identification of populations in your file (lets say that CHR_A), than you can just use subset function

ld <- read.table("plink.ld",header=T) ld$a <- (ld$BP_B-ld$BP_A) #I recomed to use meaningful names, not "a" for(family in levels(ld$CHR_A)){ pop_ld <- subset(ld, family == ld$CHR_A); pdf(paste('ld',family,'.pdf', sep = '')); plot(pop_ld$a,pop_ld$R2); dev.off() }

P.S. I am sorry, I do not know how to format code here... Just assume, that ; is a newline.

ADD REPLY
0
Entering edit mode

Thanks @kamiljaron but there are no population identifiers here. I think i need to subset the data right there in Plink while calculating the LD. So does 10 populations mean 10 different plink input .ped files?

ADD REPLY

Login before adding your answer.

Traffic: 1511 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6