Things that really don’t matter: megabase or megabasepair

Should we talk about physical distance in genetics as number of base pairs (kbp, Mbp, and so on) or bases (kb, Mb)?

I got into a discussion about this recently, and I said I’d continue the struggle on my blog. Here it is. Let me first say that I don’t think this matters at all, and if you make a big deal out of this (or whether ”data” can be singular, or any of those inconsequential matters of taste we argue about for amusement), you shouldn’t. See this blog post as an exorcism, helping me not to trouble my colleagues with my issues.

What I’m objecting to is mostly the inconsistency of talking about long stretches of nucleotides as ”kilobase” and ”megabase” but talking about short stretches as ”base pairs”. I don’t think it’s very common to call a 100 nucleotide stretch ”a 100 b sequence”; I would expect ”100 bp”. For example, if we look at Ensembl, they might describe a large region as 1 Mb, but if you zoom in a lot, they give length in bp. My impression is that this is a common practice. However, if you consistently use ”bases” and ”megabase”, more power to you.

Unless you’re writing a very specific kind of bioinformatics paper, the risk of confusion with the computer storage unit isn’t a problem. But there are some biological arguments.

A biological argument for ”base”, might be that we care about the identity of the base, not the base pairing. We note only one nucleotide down when we write a nucleic acid sequence. The base pair is a different thing: that base bound to the one on the other strand that it’s paired with, or, if the DNA or RNA is single-stranded, it’s not even paired at all.

Conversely, a biochemical argument for ”base pair” might be that in a double-stranded molecule, the base pair is the relevant informational unit. We may only write one base in our nucleotide sequence for convenience, but because of the rules of base pairing, we know the complementing pair. In this case, maybe we should use ”base” for single-stranded molecules.

If we consult two more or less trustworthy sources, The Encylopedia of Life Sciences and Wiktionary, they both seem to take this view.

eLS says:

A megabase pair, abbreviated Mbp, is a unit of length of nucleic acids, equal to one million base pairs. The term ‘megabase’ (or Mb) is commonly used inter-changeably, although strictly this would refer to a single-stranded nucleic acid.

Wiktionary says:

A length of nucleic acid containing one million nucleotides (bases if single-stranded, base pairs if double-stranded)

Please return next week for the correct pronunciation of ”loci”.

Literature

Dear, P.H. (2006). Megabase Pair (Mbp). eLS.

Reflektion om högskolepedagogik, apropå kurs

Kära dagbok,

Det finns väl inget utryck som för en att låta så mycket som en universitetslärare som ”reflektion”. Det skulle vara om ordet står i plural och stavas med x. Jag går en fortsättningskurs i högskolepedagogik just nu, och det får mig att tänka på lärande.

En positiv bieffekt är att den låter mig resa till Alnarp, ett av SLU:s andra campusområden. Det är fint där.

För det första är jag förvånad över hur mycket som sitter kvar från grundkursen jag gick som doktorand. Men för säkerhets skull, för jag hade känslan av att jag lärde mig ganska lite då, skummade jag igenom Biggs & Tang (2011), kursboken för SLU:s grundkurs i högskolepedagogik. Den håller hårt på constructive alignment av mål, aktiviteter och examination: idén att en bra kurs först definierar var studenterna ska kunna i slutet av den, och sedan övar och examinerar den just detta. Det är viktigt att precisera just vad ”kunna” betyder i sammanhanget: ska studenterna kunna förklara begrepp, ska de kunna lösa problem, eller ska de kunna bygga en maskin? I så fall bör de få öva sig på att förklara, eller bygga maskiner, och examineras på sin förmåga att förklara, eller bygga maskiner. Det är en av de där förrädiskt enkla idéerna som kan leda till dramatiska förändringar om man tar den på allvar.

Jag tänker förstås tillbaka på när jag var student. En del kurser var exemplariska i att sätta oss i basgrupper där vi fick förklara begrepp för varandra eller laborationer där vi fick lösa problem. En gästföreläsare (vars namn jag glömt; hon kom från Rättsmedicinalverket, tror jag) gav en 15 minuters introduktion till sitt ämne. Sedan lät hon oss tänka ut vad vi ville veta, formulera frågor, och sedan byggde hon resten av föreläsningen kring dem. Intressant sätt att både ge studenterna ansvar och hjälpa henne tackla vår förförståelse av ämnet. Modigt också. Annat var mindre exemplariskt; många sömniga föreläsningar blev det.

Själv har jag ofta tänkt, och sagt, ungefär ”tänk om jag kunde gå tillbaka och göra om de där första årens kurser med allt jag lärt mig om att lära mig”. I viss mån är det säkert sant att jag blev bättre på att lära mig. De sista åren som student fick jag till exempel för mig att organisera mina anteckningar som färgglada tankekartor, och repetera genom att titta på delar av dem och förklara dem för mig själv. Det var säkert inte helt tokigt. Å andra sidan, mycket av det jag lärde mig om att lära mig var nog också att strategiskt och ytligt ta in det som behövdes för att lyckas med tentan. På så sätt är det kanske till och med möjligt att lära sig sämre sämre med mer erfarenhet.

Som deltagare i en kurs om pedagogik fäster en naturligtvis extra uppmärksamhet vid hur undervisningen går till. Lever lärarna som de lär? Och vilka aktiviteter kan jag snappa upp och själv använda, helst de som inte kräver att kursen designas om i grunden utan bara annan användning av den tid och de resurser som finns? Kursen har inte gjort mig besviken så långt. Den fungerar verkligen som demonstration av sitt eget innehåll, med massor av små think–pair–share-moment, tillverkning av konceptkartor i grupper, gott om tid till diskussion, och aktiviteter som låter oss använda det vi läst in oss på i förväg, så att det får förberedelsen att kännas meningsfull.

Ett par höjdpunkter från de första kurstillfällena:

Linn Areskoug från enheten för pedagogisk utveckling samt från Uppsala universitet pratade om normkritik. Rubriken var ”könsmedveten undervisning”, men som sig bör var innehållet ifrågasättande även mot den idén. Det slog mig också under föreläsningen att föreläsningsstilen till stod del bestod av att problematisera de begrepp som just introducerats med en serie exempel och öppna frågor. När syftet var att åhörarna skulle tänka över identitet och ifrågasätta saker som känns uppenbara för oss, är det kanske en bra strategi. Jag undrar om det skulle gå att presentera svåra problem inom genetik på ett liknande sätt, med betoningen på frågorna, inte svaren.

Alexis Engström från Uppsala universitet gav oss, bland annat, bokstavligt talat en lista på aktiviteter som bryter med traditionell undervisning. Det jag tyckte lät mest spännande var ”Det saknade perspektivet”: att lämna ett pass i schemat öppet och ge studenterna i uppgift att fylla det med innehåll och aktiviteter som de tycker saknas i kursen, lite som en storskalig variant av vad gästföreläsaren från Rättsmedicinalverket gjorde.

Jag hoppas få med mig några saker jag kan göra för att vara en bättre lärare, och förstå lite bättre hur lärande fungerar. Som bonus har jag i alla fall roligt.

Litteratur

Biggs J & Tang C. (2011) Teaching for quality learning at university. 4th edition. Open University Press.

If research is learning, how should researchers learn?

I’m taking a course on university pedagogy to, hopefully, become a better teacher. While reading about students’ learning and what teachers ought to do to facilitate it, I couldn’t help thinking about researchers’ learning, and what we ought to do to give ourselves a good learning environment.

Research is, largely, learning. First, a large part of any research work is learning what is already known, not just by me in particular; it’s a direct continuation of learning that takes place in courses. While doing any research project, we learn the concepts other researchers use in this specific sub-subfield, and the relations between them. First to the extent that we can orient ourselves, and eventually to be able to make a contribution that is intelligible to others who work there. We also learn their priorities, attitudes and platitudes. (Seriously, I suspect you learn a lot about a sub-subfield by trying to make jokes about it.) We also learn to do something new: perform a laboratory procedure, a calculation, or something like that.

But more importantly, research is learning about things no-one knows yet. The idea of constructivist learning theory seems apt: We are constructing new knowledge, building on pre-existing structures. We don’t go out and read the book of nature; we take the concepts and relations of our sub-subfield of choice, and graft, modify and rearrange them into our new model of the subject.

If there is something to this, it means that old clichéd phrases like ”institution of higher learning”, scientists as ”students of X”, and so on, name a deeper analogy than it might seem. It also suggests that innovations in student learning might also be good building blocks for research group management. Should we be concept mapping with our colleagues to figure out where we disagree about the definition of ”developmental pleiotropy”? It also makes one wonder why meetings and departmental seminars often take the form of sage on the stage lectures.

Forskare vill prata om den nya gentekniken, men vem vill prata med forskare?

Forskare: ”Vi måste prata!”

Tystnad. Gryllidae-läten hörs i bakgrunden.

Andra forskare: ”Äh, håll tyst!”

För några veckor sedan kom det en debattartikel i SvD: ”Vi måste ta ställning till (eller samtala om, beroende på var på sidan man läser) den nya gentekniken” från ett gång forskare kopplade till gruppen CRISPRideas vid Pufendorfinstitutet vid Lunds universitet. Deras uppgift verkar vara att analysera hur genredigeringstekniker som CRISPR/Cas debatteras:

Genom jämförande analyser av debatten i framförallt de nordiska länderna, men även globalt, kommer vi att undersöka hur vetenskaplig kunskap, värderingar och normer har påverkat såväl experter som allmänheten och olika intressenter i deras uppfattning om de nya genredigeringsteknikernas möjligheter och risker.

Det förefaller som att debattartikeln skrevs i samband med deras möte ”Medical and agricultural perspectives on new genome editing technologies” i november, och det ser ut att ha varit ett kul möte.

Synd att debattartikeln blev tråkig. Men det är kanske inte så konstigt. Om man är en grupp forskare med olika bakgrund ens mandat är att undersöka om en debatt, så är kanske den mest naturliga åsikten att driva att en bred debatt är bra och nödvändig.

Det kan låta som en åsikt som det är svårt att invända mot, men komiskt nog finns det andra som gör det. För strax därpå tyckte Jesper Sundström och Torbjörn Fagerström att det var ”Risk för ohederlig debatt om gensaxen”.

Några kommentarer om debattartikeln. Jag förstår att de hade begränsat utrymme, men här på min blogg kan jag breda ut mig hur mycket jag vill, så vi kan titta närmare på några påståenden.

Get some of that prime jive, get some of that, get get get down

Artikeln börjar med prime-redigering, en ny variant av CRISPR/Cas-redigering som publicerades nyligen (Anzalone & al 2019). Det är tydligen på grund av den som det är extra viktigt att debattera. (Oklart vad det borde heta på svenska: ”prime” är ”början” eller ”start”, som i ”primer”, en kort RNA eller DNA-bit som hjälper en enzymatisk dna-kopiering starta.)

Idén med prime-redigering är att istället för att klippa av hela dna-strängen och klistra ihop den med en ny, så tar man bara av halva i taget, för att inte råka ha sönder för mycket, och syntetiserar det nya dna:t på plats, så att det säkert kommer i kontakt med den nyklippta dna-strängen. Det är en påhittig konstruktion som kombinerar CRISPR-systemet, stulet från bakterier, med omvänt transkriptas, stulet från virus.

Varför behövs all den här påhittigheten? För att genredigering i själva verket fungerar ganska dåligt. Att slå ut en gen i cellkultur är en relativt smal sak: klipp ett hål och låt cellen klistra ihop den bäst den kan, förmodligen blir resultatet en trasig gen. Men att byta ut en variant av en gen med en annan är svårare; majoriteten av försöken misslyckas. Det är okej i labbet, när det bara är att odla nya celler och försöka igen. Det skulle vara mindre okej i en genterapisituation där en patient ska få en ny fungerande genvariant. Prime-redigering kanske är lösningen.

Det är inte lätt att veta om prime-redigering kommer bli det nya, eller om det bara är en i serien av olika förfinade CRISPR/Cas-varianter. Som utgångspunkt för debatt ställer den inga nya frågor som inte redan ställs av tidigare CRISPR/Cas-varianter, precis som CRISPR/Cas inte leder till några nya frågor som inte redan ställdes av rekombinant DNA och fosterdiagnostik. Däremot kanske den gör att genredigering blir lättare att få att fungera, och på så sätt kan de frågorna bli mer aktuella.

Har prime-redigering botat 175 genetiska sjukdomar? Nej, ännu inte en enda.

I debattartikeln beskriver de vad prime-redigering kan göra så här. Det är föredömligt korrekt, men det är ändå värt att understryka vad det betyder:

David Liu och hans forskargrupp visar att den nya gensaxen fungerar anmärkningsvärt väl: de lyckades korrigera ett häpnadsväckande antal mutationer (175 stycken!) i olika celler från möss och människor. Bland annat har man lyckats korrigera de tidigare svåråtkomliga sjukdomsmutationerna som orsakar sickelcell-anemi och den allvarliga nervsjukdomen Tay Sachs sjukdom. Möjligheterna med tekniken är stora. Liu och hans forskargrupp hävdar att över 80 procent av alla nu kända sjukdomsalstrande mutationer kan korrigeras med prime-redigering. (min kursivering)

Det prime-redigering löser är att den gör det lättare att faktiskt redigera, inte bara ta sönder. Det har varit svårt för CRISPR/Cas i praktiken att ersätta en genvariant med en annan. Men för att göra någon form av genterapi med prime-redigering måste någon också se till att det går att rikta redigeringarna till det organ eller den vävnad där problemet finns, och visa att det är tillräckligt att den fungerande genvarianten uttrycks där för att sjukdomsförloppet ska vändas, eller åtminstone avstanna.

Är de dna-redigerade kinesiska barnen skyddade mot hiv?

Sedan pratar de om det kinesiska fallet med barn som (under falska förespeglingar) fått sitt dna-redigerat i en gen som har med hiv-infektion att göra. Den här formuleringen, däremot, är så slarvig att den inte riktig är sann:

I fjol använde en kinesisk forskare Crispr för att förändra arvsmassan hos foster så att de skyddas mot hiv. Detta ledde till en storm i forskarvärlden eftersom permanent förändring av arvsmassan, som förs vidare till kommande generationer, är förbjudet i Sverige och många andra länder.

Det finns inget snällare sätt att beskriva kunskapsläget än att den kinesiska forskaren (He Jiankui) hoppades att barnen skulle bli resistenta mot hiv. Det är inte alls säkert att de faktiskt är resistenta. Vad han gjorde var att klippa hål i genen CCR5 som kodar för en av de receptorer som hiv kan använda för att komma in i cellen. Mutationerna var inspirerade av en allel som kallas CCR5-delta32 och som ger resistens mot vissa typer av hiv genom att slå ut receptorn. Förmodligen har mutationerna som barnen bär på också den effekten att de slår ut receptorn, och i så fall är de förmodligen resistenta mot vissa typer av hiv. Dessutom är det mycket möjligt att barnen har olika mutationer i olika delar av kroppen, så de kanske har en fungerande CCR5 i vissa kroppsdelar och en trasig i andra.

Jag tycker det är fel att koncentrera sig på att det var en permanent förändring i arvsmassan, som om det var det enda problemet. Att barnen, när de växer upp, kommer behöva oroa sig för om eventuella problem går i arv är illa, men det är bara ett i en rad av risker de utsätts för. Det är frågan om en oprövad behandling, med helt okända biverkningar, dåligt genomförd (fråga vem som helst som jobbat med CRISPR/Cas själv vad de tycker om Jiankuis kvalitetskontroll av redigeringar), utförd på ofödda barn vars föräldrar nästan säkert fått bristfällig information … Att lägga fram det som om det var frågan om en behandling som faktiskt gjorde nytta, men som tyvärr är ärftlig är att missa poängen.

Finns det någon risk att Sverige hamnar på efterkälken pga av etiska kval?

Författarna oroar sig, efter det kinesiska fallet, att genredigering ska förbjudas och Sverige på något sätt hamna på efterkälken. Det hade varit hjälpsamt om de givit något exempel på vem i Sverige som tycker detta?

Eftersom det finns en stark opinion mot GMO och det är svårt att få tillstånd för odling av GMO-baserade grödor i Europa har den nya regleringen kraftigt bromsat forskningen inom området. Det hindrar utvecklingen av hållbara grödor. Skulle samma sak drabba Crispr-tekniker för medicinsk tillämpning kommer den vetenskapliga utvecklingen att stanna av och våra möjligheter att bota allvarliga genetiska sjukdomar kommer att reduceras.

Genetiska förändringar i människor som går i arv är inte tillåtet i Sverige, oavsett teknik (Lag (2006:351) om genetisk integritet m.m.). Som synes är lagen från 2006, så inget nytt under solen.

Debatt om debatten om debatten

Till sist, det här är kanske långsökt, men debatten om debatten fick mig att tänka på detta fina referat av svensk lärd debatt från Tage Danielssons Grallimatik (1966).

Professor Gunnar Biörck:
Jag gör mig härmed till tolk för en tigande opinion som ogillar normlöshet och trolöshet och den junta av kulturradikaler som skriver om sånt.

Fil. dr Olof Lagercrantz:
Det är farligt med okunniga professorer.

Fil. lic. Johan Asplund:
Naturvetenskaparna borde någon gång säga något.

Naturvetenskaparna borde någon gång säga något. Det tycker undertecknarna också:

I vårt tvärvetenskapliga projekt vid Pufendorfinstitutet på Lunds universitet – CRISPRideas – bryts naturvetenskapliga och medicinska synsätt mot etiska, filosofiska, juridiska och humanistiska synsätt. Det har vidgat våra vyer och stärkt oss i vår uppfattning att det vetenskapliga samfundet kan och skall vara en partner i svåra men viktiga samtal. Björnen har lämnat sitt ide – den är inte farlig, bara man är varlig! (Slutklämmen på slutrepliken, ”Att inte debattera är inte ett alternativ”.)

Litteratur

Antonio Regalado. (2019) China’s CRISPR: babies. MIT Technology Review.

Andrew V Anzalone et al. (2019) Search-and-replace genome editing without double-strand breaks or donor DNA. Nature.

Two distinguishing traits of science are that there are errors all the time and that almost no-one can reproduce anything

I got annoyed and tweeted:

”If you can’t reproduce a result, it isn’t science” … so we’re at that stage now, when we write things that sound righteous but are nonsense.

Hashtag subtweet, I guess. But it doesn’t matter who first wrote the sentence I was complaining about; they won’t care what I think, and I’m not out to debate them. I only think the quoted sentence makes sense if you take ”science” to mean ”the truth”. The relationship between science and reproducibility is messier than that.

The first clause could mean a few different things:

You have previously produced a result, but now, you can’t reproduce it when you try …

Then you might have done something wrong the first time, or the second time. This is an everyday occurrence of any type of research, that probably happens to every postdoc every week. Not even purely theoretical results are safe. If the simulation is stochastic, one might have been interpreting noise. If there is an analytical result, one might have made an odd number of sign errors. In fact, it is a distinguishing trait of science that when we try to learn new things, there are errors all the time.

If that previous result is something that has been published, circulated to peers, and interpreted as if it was a useful finding, then that is unfortunate. The hypothetical you should probably make some effort to figure out why, and communicate that to peers. But it seems like a bad idea to suggest that because there was an error, you’re not doing science.

You personally can’t reproduce a results because you don’t have the expertise or resources …

Science takes a lot of skill and a lot of specialised technical stuff. I probably can’t reproduce even a simple organic chemistry experiment. In fact, it is a distinguishing trait of science that almost no-one can reproduce any of it, because it takes both expertise and special equipment.

No-one can ever reproduce a certain result even in principle …

It might still be science. The 1918 influenza epidemic will by the nature of time never happen again. Still, there is science about it.

You can’t reproduce someone else’s results when you try with a reasonably similar setup …

Of course, this is what the original authors of the sentence meant. When this turns out to be a common occurrence, as people systematically try to reproduce findings, there is clearly something wrong with the research methods scientists use: The original report may be the outcome of a meandering process of researcher degrees of freedom that produced a striking result that is unlikely to happen when the procedure is repeated, even with high fidelity. However, I would say that we’re dealing with bad science, rather than non-science. Reproducibility is not a demarcation criterion.

(Note: Some people reserve ”reproducibility” for the computational reproducibility of re-running someone’s analysis code and getting the same results. This was not the case with the sentence quoted above.)

Self-indulgent meta-post of the year

Dear diary,

Time for a recap of the On unicorns and genes blogging year. During 2019, your friendly neighbourhood genetics blog mostly kept to its schedule of four posts per month with some blog vacation in summer and in December.

This added up to a total of 43 posts (including this one), only one of them in Swedish (Gener påverkar ditt och datt, reflecting on how genome-wide association is often reported as if it was something else), and three of them posts about three first-author papers that came out in 2019:

Paper: ”Integrating selection mapping with genetic mapping and functional genomics”
Paper: ”Sequence variation, evolutionary constraint, and selection at the CD163 gene in pigs”
Paper: ”Removal of alleles by genome editing (RAGE) against deleterious load”

Now, let’s pick one post per month to represent the blogging year of 2019:

January: Showing a difference in means between two groups. This is one of those hard easy problems: Imagine we have an experiment comparing the means of two groups; how do we show both the data and the estimate of the difference, with uncertainty, in the same plot? In this little post, I try a more and a less radical version. I’ve since used the less radical version a couple of times.

February: ”We have reached peak gene and passed it”. Comment on an opinion piece that argued that revisions to the gene concept have important implications for modern genetics, and people need to be told. I agreed about a lot of the criticisms, but thought they should have nothing to do with gene concepts.

March: Journal club of one: ”Biological relevance of computationally predicted pathogenicity of noncoding variants”. Journal club post about a then recent paper about how variant effect prediction, especially for noncoding variants, is really hard, and not easy to evaluate fairly either.

April: Greek in biology. Comment on a stimulating essay about multilingualism in biology.

May: What single step does with relationship. A simulation and a couple of heatmaps to try to understand how the single step method describes relationship between individuals by blending genomic and pedigree relatedness.

June: Simulating genetic data with R: an example with deleterious variants Post from my talk at the Edinburgh R users group.

July: Using R: Correlation heatmap with ggplot2. I thought I’d update one of my most viewed posts, which was becoming too embarrassingly outdated. Unfortunately, I also broke its old address, and that wasn’t so smart.

August: Blog vacation.

September: Genes do not form networks. They just don’t. It seems DiFrisco & Jaeger (2019) don’t think so either. I would like to return to this later.

October: Using R: Animal model with simulated data. What it says on the tin: example of simulating a pedigree with AlphaSimR and fitting an animal model.

November: Using R: from gather to pivot. Introduction to, and celebration of, one of the best changes in the tidyverse.

December: Interpreting genome scans, with wisdom. Opinions about genome-wide association, precipitated by Eric Faumann’s Twitter account.

This is also the year I moved from the Roslin in Scotland to Uppsala, Sweden, for the second phase of the mobility grant-supported postdoc. Here are many of my worldly possessions, being hauled:

What will happen on here in 2020? The ambition is to keep a reasonably regular schedule, post about papers as they come out, and maybe write some more in Swedish.

I’ll also have a blog anniversary. The blog, admittedly rather different the first years, started in earnest in June 2010. I’m not sure how to celebrate, but I feel like I should follow up on one of the first posts (about linkage mapping and the probably spurious ”Gay Gene” at Xq28), now that Ganna et al. (2019) put genetic mapping of sexual orientation is in the news again.

Using R: Animal model with hglm and Stan (with Cholesky trick)

A few weeks ago I posted about fitting the quantitative genetic animal model with MCMCglmm and R-INLA. Since then, I listened to a talk by Lars Rönnegård, one of the creators of the hglm package, and this paper was published in GSE about animal models in Stan.

hglm

The hglm package fits hierarchical generalised linear models. That includes the animal model with pedigree or genomic relatedness. Hierarchical generalised linear models also allow you to model the dispersion of random effects, which lets you do tricks like variance QTL mapping (Rönnegård & Valdar 2011), breeding values for variances (Rönnegård et al. 2010) or genomic prediction models with predictors of marker variance (Mouresan, Selle & Rönnegård 2019). But let’s not get ahead of ourselves. How do we fit an animal model?

Here is the matrix formulation of the animal model that we skim through in every paper. It’s in this post because we will use the design matrix interface to hglm, which needs us to give it these matrices (this is not a paper, so we’re not legally obliged to include it):

\mathbf{y} = \mu + \mathbf{X} \mathbf{b} + \mathbf{Z} \mathbf{a} + \mathbf{e}

The terms are the the trait value, intercept, fixed coefficients and their design matrix, genetic coefficients and their design matrix, and the residual. The design matrix Z will contain one row and column for each individual, with a 1 to indicate its position in the phenotype table and pedigree and the rest zeros. If we sort our files, it’s an identity matrix.

The trick with the genetic coefficients is that they’re correlated, with a specific known correlation structure that we know from the pedigree (or in genomic models, from markers). It turns out (Lee, Nelder & Pawitan 2017, chapter 8) that you can change the Z matrix around so that it lets you fit the model with an identity covariance matrix, while still accounting for the correlations between relatives. You replace the random effects for relatedness with some transformed random effects that capture the same structure. One way to do this is with Cholesky decomposition.

\mathbf{Z_{fudged}} = \mathbf{Z_0} \mathbf{L}

As an example of what the Cholesky decomposition does, here is slice of the additive relationship matrix of 100 simulated individuals (the last generation of one replicate of these simulations) and the resulting matrix from Cholesky decomposition.

So instead of

\mathbf{a} \sim N(0, \mathbf{A} \sigma)

We can fit

\mathbf{a_{fudged}} \sim N(0, \mathbf{I} \sigma)

This lets us fit the animal model with hglm, by putting in a modified Z matrix.

Assuming we have data frames with a pedigree and a phenotype (like, again, from these simulations):

library(AGHmatrix)
library(hglm)

A  <- Amatrix(ped)

Z0  <- diag(1000)
L <- t(chol(A))
Z  <- Z0 %*% L
X <- model.matrix(~1, pheno)

model <- hglm(y = pheno$pheno,
              X = X,
              Z = Z,
              conv = 1e-8)

est_h2  <- model$varRanef / (model$varRanef + model$varFix)

(I found the recommendation to decrease the convergence criterion from the default for animal models in a YouTube video by Xia Chen.)

Stan

When we turn to Stan, we will meet the Cholesky trick again. Stan is a software for Markov Chain Monte Carlo, built to fit hierarchical linear models, and related high-dimensional models, more effectively than other sampling strategies (like Gibbs). rstan is a helpful package for running Stan from within R.

Nishio & Arakawa (2019) recently published a Stan script to fit an animal model, comparing Stan to a Gibbs sampler (and a related MCMC sampler that they also didn’t publish the code for). If we look into their Stan model code, they also do a Cholesky decomposition to be able to use an identity matrix for the variance.

First, they decompose the additive relationship matrix that the program takes in:

transformed data{
  matrix[K,K] LA;
  LA = cholesky_decompose(A);
}

And then, they express the model like this:

vector[N] mu;
vector[K] a;
a_decompose ~ normal(0, 1);
a = sigma_G * (LA * a_decompose);
mu = X * b + Z * a;
Y ~ normal(mu, sigma_R);

We can add this line to the generated quantities block of the Stan program to get heritability estimates directly:

real h2;
h2 = sigma_U / (sigma_U + sigma_E)

Here, we’ve saved their model to a stan file, and now we can run it from R:

pheno$scaled_pheno <- as.vector(scale(pheno$pheno))

model_stan <- stan(file = "nishio_arakawa.stan",
                   data = list(Y = pheno$scaled_pheno,
                               X = X,
                               A = A,
                               Z = Z0,
                               J = 1,
                               K = 1000,
                               N = 1000))

est_h2_stan <- summary(model_stan, pars = "h2")$summary

Important note that I always forget: It's important to scale your traits before you run this model. If not, the priors might be all wrong.

The last line pulls out the summary for the heritability parameter (that we added above). This gives us an estimate and an interval.

The paper also contains this entertaining passage about performance, which reads as if it was a response to a comment, actual or anticipated:

R language is highly extensible and provides a myriad of statistical and graphical techniques. However, R language has poor computation time compared to Fortran, which is especially well suited to numeric computation and scientific computing. In the present study, we developed the programs for GS and HMC in R but did not examine computation time; instead, we focused on examining the performance of estimating genetic parameters and breeding values.

Yes, two of their samplers (Gibbs and HMC) were written in R, but the one they end up advocating (and the one used above), is in Stan. Stan code gets translated into C++ and then compiled to machine code.

Stan with brms

If rstan lets us run Stan code from R and examine the output, brms lets us write down models in relatively straightforward R syntax. It’s like the MCMCglmm of the Stan world. We can fit an animal model with brms too, by directly plugging in the relationship matrix:

model_brms <- brm(scaled_pheno ~ 1 + (1|animal),
                  data = pheno,
                  family = gaussian(),
                  cov_ranef = list(animal = A),
                  chains = 4,
                  cores = 1,
                  iter = 2000)

Then, we can pull out the posterior samples for the variability, here expressed as standard errors, compute the heritability and then get the estimates (and interval, if we want):

posterior_brms <- posterior_samples(model_brms,
                                    pars = c("sd_animal", "sigma"))

h2_brms  <- posterior_brms[,1]^2 /
    (posterior_brms[,1]^2 + posterior_brms[,2]^2)

est_h2_brms <- mean(h2_brms)

(Code is on GitHub: both for the graphs above, and the models.)

Interpreting genome scans, with wisdom

Eric Fauman is a scientist at Pfizer who also tweets out interpretations of genome-wide association scans.

Background: There is a GWASbot twitter account which posts Manhattan plots with links for various traits from the UK Biobank. The bot was made by the Genetic Epidemiology lab at the Finnish Institute for Molecular Medicine and Harvard. The source of the results is these genome scans (probably; it’s little bit opaque); the bot also links to heritability and genetic correlation databases. There is also an EnrichrBot that replies with enrichment of chromatin marks (Chen et al. 2013). Fauman’s comments on some of the genome scans on his Twitter account.

Here are a couple of recent ones:

And here is his list of these threads as a Google Document.

This makes me thing of three things, two good, and one bad.

1. The ephemeral nature of genome scans

Isn’t it great that we’re now at a stage where a genome scan can be something to be tweeted or put en masse in a database, instead of published one paper per scan with lots of boilerplate. The researchers behind the genome scans say as much in their 2017 blog post on the first release:

To further enhance the value of this resource, we have performed a basic association test on ~337,000 unrelated individuals of British ancestry for over 2,000 of the available phenotypes. We’re making these results available for browsing through several portals, including the Global Biobank Engine where they will appear soon. They are also available for download here.

We have decided not to write a scientific article for publication based on these analyses. Rather, we have described the data processing in a detailed blog post linked to the underlying code repositories. The decision to eschew scientific publication for the basic association analysis is rooted in our view that we will continue to work on and analyze these data and, as a result, writing a paper would not reflect the current state of the scientific work we are performing. Our goal here is to make these results available as quickly as possible, for any geneticist, biologist or curious citizen to explore. This is not to suggest that we will not write any papers on these data, but rather only write papers for those activities that involve novel method development or more complex analytic approaches. A univariate genome-wide association analysis is now a relatively well-established activity, and while the scale of this is a bit grander than before, that in and of itself is a relatively perfunctory activity. [emphasis mine] Simply put, let the data be free.

That being said, when starting to write this post, first I missed a paper. It was pretty frustrating to find a detailed description of the methods: after circling back and forth between the different pages that link to each other, I landed on the original methods post, which is informative, and written in a light conversational style. On the internet, one would fear that this links may rot and die eventually, and a paper would probably (but not necessarily …) be longer-lasting.

2. Everything is a genome scan, if you’re brave enough

Another thing that the GWAS bot drives home is that you can map anything that you can measure. The results are not always straightforward. On the other hand, even if the trait in question seems a bit silly, the results are not necessarily nonsense either.

There is a risk, for geneticists and non-geneticists alike, to reify traits based on their genetic parameters. If we can measure the heritability coefficient of something, and localise it in the genome with a genome-wide association study, it better be a real and important thing, right? No. The truth is that geneticists choose traits to measure the same way all researchers choose things to measure. Sometimes for great reasons with serious validation and considerations about usefulness. Sometimes just because. The GWAS bot also helpfully links to the UK Biobank website that describes the traits.

Look at that bread intake genome scan above. Here, ”bread intake” is the self-reported number of slices of bread eaten per week, as entered by participants on a touch screen questionnaire at a UK Biobank assessment centre. I think we can be sure that this number doesn’t reveal any particularly deep truth about bread and its significance to humanity. It’s a limited, noisy, context-bound number measured, I bet, because once you ask a battery of lifestyle questions, you’ll ask about bread too. Still, the strongest association is at a region that contains olfactory receptor genes and also shows up two other scans about food (fruit and ice cream). The bread intake scan hits upon a nugget of genetic knowledge about human food preference. A small, local truth, but still.

Now substitute bread intake for some more socially relevant trait, also imperfectly measured.

3. Lost, like tweets in rain

Genome scan interpretation is just that: interpretation. It means pulling together quantitative data, a knowledge of biology, previous literature, and writing an unstructured text, such as a Discussion section or a Twitter thread. This makes them harder to organise, store and build on than the genome scans themselves. Sure, Fauman’s Twitter threads are linked from the above Google Document, and our Discussion sections are available from the library. But they’re spread out in different places, they mix (as they should) evidence with evaluation and speculation, and it’s not like we have a structured vocabulary for describing genetic mechanisms of quantitative trait loci, and the levels of evidence for them. Maybe we could, with genome-wide association study ontologies and wikis.

You’re not funny, but even if you were

Here is a kind of humour that is all too common in scientific communication; I’ll just show you the caricature, and I think you’ll recognize the shape of it:

Some slogan about how a married man is a slave or a prisoner kneeling and holding a credit card. Some joke where the denouement relies on: the perception that blondes are dumb, male preference for breast size, perceived associations between promiscuity and nationality, or anything involving genital size. Pretty much any one-panel cartoon taken from the Internet.

Should you find any of this in your own talk, here is a message to you: That may be funny to you; that isn’t the problem. To a fair number of the people who are listening, it’s likely to be trite, sad and annoying.

Humour totally has a place in academic speech and writing—probably more than one place. There is the laughter that is there to relieve tension. That is okay sometimes. There are jokes that are obviously put-downs. Those are probably only a good idea in private company, or in public forums where the object of derision is powerful enough that you’re not punching down, but powerless enough to not punch you back. Say, the ever-revered and long dead founder of your field—they may deserve a potshot at their bad manners and despicable views on eugenics.

Then there is that elusive ‘sudden perception of the incongruity between a concept and the real objects which have been thought through it in some relation’ (Schopenhauer, quoted in Stanford Encyclopedia of Philosophy). When humour is used right, a serious lecturer talking about serious issues has all kinds of opportunities to amuse the listener with incongruities between the expectations and what they really are like. So please don’t reveal yourself to be predictably trite.

Using R: Installing GenABEL and RepeatABEL

GenABEL is an R package for performing genome-wide association with linear mixed models and a genomic relationship matrix. RepeatABEL is a package for such genome-wide association studies that also need repeated measures.

Unfortunately, since 2018, GenABEL is not available on CRAN anymore, because of failed checks that were not fixed. (Checks are archived on CRAN, but this means very little to me.) As a consequence, RepeatABEL is also missing.

Fair enough, the GenABEL creators probably aren’t paid to maintain old software. It is a bit tragic, however, to think that in 2016, GenABEL was supposed to be the core of a community project to develop a suite of genomic analysis packages, two years before it was taken of CRAN:

The original publication of the GenABEL package for statistical analysis of genotype data has led to the evolution of a community which we now call the GenABEL project, which brings together scientists, software developers and end users with the central goal of making statistical genomics work by openly developing and subsequently implementing statistical models into user-friendly software.

The project has benefited from an open development model, facilitating communication and code sharing between the parties involved. The use of a free software licence for the tools in the GenABEL suite promotes quick uptake and widespread dissemination of new methodologies and tools. Moreover, public access to the source code is an important ingredient for active participation by people from outside the core development team and is paramount for reproducible research. Feedback from end users is actively encouraged through a web forum, which steadily grows into a knowledge base with a multitude of answered questions. Furthermore, our open development process has resulted in transparent development of methods and software, including public code review, a large fraction of bugs being submitted by members of the community, and quick incorporation of bug fixes.

I have no special insight about the circumstances here, but obviously the situation is far from ideal. You can still use the packages, though, with a little more effort to install. Who knows how long that will be the case, though. In a complex web of dependencies like the R package ecosystem, an unmaintained package probably won’t last.

GenABEL can probably be replaced by something like GEMMA. It does mixed models for GWAS, and while it isn’t an R package, it’s probably about as convenient. However, I don’t know of a good alternative to RepeatABEL.

These are the steps to install GenABEL and RepeatABEL from archives:

  1. We go to the CRAN archive and get the tarballs for GenABEL, GenABEL.data which it needs, and RepeatABEL.
    curl -O https://cran.r-project.org/src/contrib/Archive/GenABEL/GenABEL_1.8-0.tar.gz
    curl -O https://cran.r-project.org/src/contrib/Archive/GenABEL.data/GenABEL.data_1.0.0.tar.gz
    curl -O https://cran.r-project.org/src/contrib/Archive/RepeatABEL/RepeatABEL_1.1.tar.gz
    

    We don’t need to unpack them.

  2. Install GenABEL.data and GenABEL from a local source. Inside R, we can use install.packages, using the files we’ve just downloaded instead of the online repository.
    install.packages(c("GenABEL.data_1.0.0.tar.gz", "GenABEL_1.8-0.tar.gz"), repos = NULL)
    
  3. To install RepeatABEL, we first need hglm, which we can get from CRAN. After that has finished, we install RepeatABEL, again from local source:
    install.packages("hglm")
    install.packages("RepeatABEL_1.1.tar.gz", repos = NULL)
    

This worked on R version 3.6.1 running on Ubuntu 16.04, and also on Mac OS X.

Literature

Karssen, Lennart C., Cornelia M. van Duijn, and Yurii S. Aulchenko. ”The GenABEL Project for statistical genomics.” F1000Research 5 (2016).