Entering edit mode
8.7 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)
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
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
functionld <- 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?