R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

ABO-veriryhmän alleelifrekvenssit

Hardy-Weinbergin laki
Hardy-Weinbergin laki (HW) kuvaa eri genotyyppifrekvenssien ja alleelifrekvenssien suhdetta vapaasti sekoittuvassa (tasapaino)populaatiossa. Laki yksinkertaistaa eri genotyyppien frekvenssit kaavaksi 1 = (p+q)^2 = p^2+2pq+q^2, jossa p on tietyn lokuksen dominoivan alleelin frekvenssi ja q resessiivisen alleelin frekvenssi populaatiossa. Kaava lienee useimmille tuttu jo lukiosta. Esitetty kaava ei sellaisenaan kuvaa tilannetta kolmen alleelin lokuksessa, jollainen ABO-veriryhmän tapaus on. ABO-lokuksessa on kolme alleelia, A, B ja o, joista A ja B ovat kodominantteja ja o on resessiivinen (oikeasti alleeleja tunnetaan kyllä useampia). Tällöin tasapainotilannetta kuvaa kaava 1 = (p+q+r)^2 = p^2+q^2+r^2+2pq+2pr+2qr.

ABO-veriryhmä ja HW
Veriryhmien ja alleelikombinaatioiden vastaavuus on seuraavanlainen:

 o = oo
 A = AA tai Ao
 B = BB tai Bo
AB = AB

Käyttäen yllä johdetun tasapainotilanteen kaavaa voidaan edelleen merkitä:

 o = oo        = r^2
 A = AA tai Ao = p^2+2pq
 B = BB tai Bo = q^2+2qr
AB = AB = pq   = 2pq

Esimerkiksi Suomessa veriryhmäjakauma on seuraava:

 o = 0.31
 A = 0.44
 B = 0.17
AB = 0.08

Tästä voidaan helposti laskea o-alleelin frekvenssi:

o = 0.31 = r^2
 <=>   r = 0.557

A-alleelin frekvenssi saadaan kaavasta:

p^2+2pr = 0.44

Sijoitetaan kaavaan edellinen tulos, ja saadaan:

            p^2 + 2p * 0.557 = 0.44
=> p2 + 2 * p * 0.557 - 0.44 = 0

Tämä on toisen asteen yhtälö p:n suhteen, joka voidaan ratkaista R:ssä vaikkapa graafisesti:

 
f<-function(p) {
   p^2 + 2 * 0.557 * p - 0.44
}
 
s<-seq(-2,2,0.001)
r<-c()
for(i in 1:length(s)) {
   print(i)
   r<-c(r, f(s[i]))
}
 
plot(x=s, y=r)
d<-data.frame(cbind(r,s))
d[d$r<=0.001 & d$r>=-0.001,]
#              r      s
# 578  -0.000293 -1.423
# 2310 -0.000293  0.309

Eli A:n frekvenssi populaatiossa on 0.309. Negatiivinen tulos jätetään huomiotta, koska frekvenssi ei voi olla negatiivinen.

Vastaavalla tavalla voitaisiin laskea myös B-alleelin frekvenssi, mutta koska frekvenssien pitää summautua ykköseen, on B:n frekvenssi siis 1 – 0.557 – 0.309 = 0.134.

Laskentavaiheista voidaan koostaa yksi funktio, joka palauttaa kaikki tulokset. Seuraavassa funktiossa on käytetty graafisen ratkaisun sijaan funktiota uniroot(), joka ratkaisee yhtälön yhden muuttujan suhteen:

hw3<-function(a, b, o) {
   r<-sqrt(o)
   f<-function(x) {
      x^2+2*x*r-a
   }
   p<-uniroot(f, c(-0.5,1))$root
   q<-1-(p+r)
   res<-c(p,q,r)
   return(res)
}

Jos funktiolle hw3() syöttää A-fenotyypin, B-fenotyypin ja o-fenotyypin frekvenssit populaatiossa, funktio palauttaa alleelien frekvenssit:

hw3(0.44,0.17,0.31)
# [1] 0.3092629 0.1339606 0.5567764

Sehän kävi näppärästi, yllätyin oikein itsekin! Mutta mitäs tällä nyt sitten tehdään? No kuvataan vaikkapa alleelifrekvenssejä eri maissa. SAS on blogissaan kuvannut fenotyyppifrekvenssien vaihtelua eri maiden välillä (1 ja 2), mutta voisi olla kiinnostavaa nähdä, miten alleellifrekvenssit muuttuvat geografisesti.

ABO-veriryhmän alleelifrekvenssit kartalla
Haetaanpa ensin ABO-veriryhmän yleisyystiedot. Ne löytyvät kahdelta sivulta html-taulukkona, joten data voidaan lukea suoraan XLM-paketin funtkiolla readHTMLTable():

library(XML)
tab1<-readHTMLTable("http://www.rhesusnegative.net/themission/bloodtypefrequencies/")
tab2<-readHTMLTable("http://en.wikipedia.org/wiki/Blood_type_distribution_by_country")

Taulukot vaativat hieman kaunistelemista, sillä kunkin veriryhmän frekvenssin perässä on prosenttimerkki, ja Rh-positiiviset ja -negatiiviset ABO-frekvenssit on eroteltu. Korjataan kyseiset ongelmat:

f<-function(x) {
   as.numeric(gsub(",", ".", gsub("%", "", x)))
}
 
tab1.1<-data.frame(tab1[[1]]$Country, apply(tab1[[1]][,3:10], 2, f))
tab1.2<-data.frame(country=tab1[[1]]$Country, O=rowSums(tab1.1[,c(2,6)]),
                   A=rowSums(tab1.1[,c(3,7)]), B=rowSums(tab1.1[,c(4,8)]),
                   AB=rowSums(tab1.1[,c(5,9)]))
tab1.3<-tab1.2[-c(66:68),]
 
tab2.1<-data.frame(tab2[[2]]$Country, apply(tab2[[2]][,3:10], 2, f))
tab2.2<-data.frame(country=tab2[[2]]$Country, O=rowSums(tab2.1[,c(2,6)]),
                   A=rowSums(tab2.1[,c(3,7)]), B=rowSums(tab2.1[,c(4,8)]),
                   AB=rowSums(tab2.1[,c(5,9)]))
tab2.3<-tab2.2[-c(21,31,32),]
tab2.3$country<-gsub("]", "", gsub("\\[", "", gsub("[[:digit:]]{1,2}", "", tab2.3$country)))

Yhdistetäänpä vielä taulukot yhdeksi isoksi taulukoksi ja poistetaan saman rivin useammat kopiot:

tab<-rbind(tab1.3, tab2.3)
tab<-tab[order(tab$country),]
tab<-tab[!duplicated(tab$country),]

Käytännössä molemmissa taulukoissa näytti olevan sama sisältö. Taulukon perusteella voidaan laskea kunkin maan alleelifrekvenssit käyttäen aiempaa funktiota:

alfreq<-data.frame(matrix(ncol=4, nrow=nrow(tab), data=NA))
alfreq[,1]<-tab$country
colnames(alfreq)<-c("country", "A", "B", "O")
for(i in 1:nrow(tab)) {
   alfreq[i,2:4]<-hw3(a=tab$A[i]/100, b=tab$B[i]/100, o=tab$O[i]/100)
}

Laskettujen alleelifrekvenssien perusteella on helpointa piirtää karttoja paketilla rworldmap:

library(rworldmap)
# Yhdistää datan karttapohjaan maiden nimien perusteella
kartta<-joinCountryData2Map(alfreq, joinCode="NAME", nameJoinColumn="country")
# Piirretään kolme karttaa, yksi kutakin alleelia kohden
png(file="ABO.png", width=530, height=1590)
par(mfrow=c(3,1))
mapCountryData(kartta, nameColumnToPlot="A", colourPalette=brewer.pal(9, "Reds"), numCats=9)
mapCountryData(kartta, nameColumnToPlot="B", colourPalette=brewer.pal(9, "Reds"), numCats=9)
mapCountryData(kartta, nameColumnToPlot="O", colourPalette=brewer.pal(9, "Reds"), numCats=9)
dev.off()

Tuloksena on alla oleva kuva, jonka tulkinta on tietysti pitkälti samanlainen kuin aiemmin mainitun SAS:n blog-postauksen; A-alleelin frekvenssi on Suomessa samalla tasolla muiden Pohjoismaiden kanssa, mutta B-alleeli on huomattavasti runsaampi. Tämä kielinee osaltaan Suomen asutushistoriasta.





Piirretäänpä vielä NMDS-kartta alleelifrekvenssien perusteella:

library(MASS)
library(RColorBrewer)
dat<-kartta@data[!is.na(kartta@data$A),]
mds<-isoMDS(dist(dat[,51:53]+rnorm(65*3, 0.00001, 0.0001)))
cols<-brewer.pal(6, "Set1")
names(cols)<-names(table(dat$REGION))[table(dat$REGION)!=0]
plot(mds$points, col="grey90", pch=15, cex=0.5)
text(x=mds$points[,1], y=mds$points[,2], labels=dat$ISO_A2, cex=0.75, col=cols[as.character(dat$REGION)], font=2)
legend(x="topleft", fill=cols, legend=names(cols), ncol=1, bty="n")
#identify(mds$points[,1], mds$points[,2], labels=dat$ISO_A2, cex=0.75)

NMDS-kartan perusteella voi yrittää hahmottaa, mitkä maat ovat keskenään kaikkein samanlaisimpia, sillä ne sijaitsevat kartassa lähimpänä toisiaan. Kartta näyttää seuraavalta, ja suurin osa Eurooppaa (oikea alareuna) ja Aasiaa (oikea yläreuna) erottuvat omiksi alueikseen, mutta muut menevät suurelta osn päällekkäin. Kuva on varsin samanlainen kuin SAS:n blogissa PCA:lla tuotettu kaavio.