This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Population-wise decay of linkage disequilibrium (LD)

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

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

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

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.

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?

0 answers

No answers yet.

Log in to answer this question.