Hi,
I have this code:
library(SNPRelate)
# get PLINK output
plink.genome <- read.table("plink.genome", header=TRUE)
> head(plink.genome)
FID1 IID1 FID2 IID2 RT EZ Z0 Z1 Z2 PI_HAT PHE DST
1 fam1054 G1054 fam1054 G700 OT 0 0.0045 0.9938 0.0017 0.4986 -1 0.839150
2 fam1054 G1054 fam1054 G701 OT 0 0.0000 1.0000 0.0000 0.5000 -1 0.838381
3 fam1079 G1079 fam2484 G2484 UN NA 0.0000 0.0007 0.9993 0.9997 -1 0.999889
4 fam1245 G1237 fam1245 G1245 OT 0 0.0036 0.9964 0.0000 0.4982 -1 0.838770
5 fam1245 G1241 fam1245 G1245 OT 0 0.0042 0.9854 0.0104 0.5031 -1 0.840569
6 fam0176 G174 fam0176 G176 OT 0 0.0000 1.0000 0.0000 0.5000 -1 0.837799
> head(plink.genome$IID1)
[1] G1054 G1054 G1079 G1237 G1241 G174
33 Levels: G1054 G1079 G1237 G1241 G174 G175 G177 G178 G1818 G2007 ... G578
snpgdsBED2GDS("output4.bed", "output4.fam","output4.bim", "HapMap.gds")
genofile <- snpgdsOpen("HapMap.gds")
# get SNPRelate output
ibd <- snpgdsIBDMoM(genofile, remove.monosnp=FALSE, kinship=TRUE)
head(ibdlist <- snpgdsIBDSelection(ibd))
ID1 ID2 k0 k1 kinship
1 G1000 G1001 1.0000000 0.00000000 0.000000000
2 G1000 G1003 0.9938901 0.00000000 0.003054932
3 G1000 G1005 1.0000000 0.00000000 0.000000000
4 G1000 G1009 1.0000000 0.00000000 0.000000000
# adjust for the orders of sample pair
pair.samp <- paste(ibdlist$ID1, ibdlist$ID2, sep=" ")
head(pair.samp)
[1] "G1000 G1001" "G1000 G1003" "G1000 G1005" "G1000 G1009" "G1000 G1052"
[6] "G1000 G1054"
plink.genome <- plink.genome[match(
paste(plink.genome$IID1, plink.genome$IID2, sep=" "), pair.samp), ]
> head(plink.genome)
FID1 IID1 FID2 IID2 RT EZ Z0 Z1 Z2 PI_HAT PHE DST PPC RATIO IBS0 IBS1
NA <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
NA.1 <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
NA.2 <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
NA.3 <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
NA.4 <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
NA.5 <NA> <NA> <NA> <NA> <NA> NA NA NA NA NA NA NA NA NA NA NA
So nothing is matching here. Can you please advise,
Thanks
Ana
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
