Jonas Giek

(1)

TVE 16 027 maj

Examensarbete 15 hp Juni 2016

Simulerad effektivisering av

genotypdataanalys genom poolade data

Simon Strömstedt Hallberg

Jonas Giek

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Simulated optimization of genotype data analysis by pooling data

Simon Strömstedt Hallberg, Jonas Giek

Målet med projektet är att undersöka om det går att effektivisera hur man

undersöker människors gener. Detta görs genom att skapa ett program i Java.

Resultatet är ett program som sorterar genotypdata från 1000 Genomes Project och utvärderar nyttan av att undersöka genotyper från flera individer

samtidigt.

Ämnesgranskare: Mia Sterby Handledare: Carl Nettelblad

(3)

Innehåll

1 Introduktion 5

1.1 Bakgrund . . . 5

1.2 Motivering och målbeskrivning . . . 5

2 Metod 6 2.1 Utförande . . . 6

2.1.1 ReadVCF . . . 6

2.1.2 RandomGroups . . . 7

2.1.3 Poolning . . . 7

2.2 Verktyg . . . 9

2.2.1 MATLAB . . . 9

2.2.2 Java . . . 9

2.2.3 HTSJDK . . . 9

2.2.4 Eclipse . . . 9

2.2.5 GitHub . . . 10

2.2.6 Maven . . . 10

3 Resultat 11 3.1 Andel entydigt bestämda individer . . . 11

3.2 Minor Allele Frequency (MAF) . . . 13

4 Diskussion 15 4.1 Utvärdering . . . 15

4.2 Felkällor . . . 16

5 Slutsatser 16 5.1 Upplevelse av verktyg . . . 16

5.1.1 Överblick . . . 16

5.1.2 GitHub . . . 16

5.1.3 Maven . . . 17

5.1.4 HTSJDK . . . 17

5.2 Utveckling . . . 17

6 Referenser 18

Appendices 19

(4)

A ReadVCF 19

B RandomGroups 21

(5)

Populärvetenskaplig sammanfattning av projektet

Grundtanken med detta projekt är att det borde gå att effektivisera hur man analyserar egenskaper hos gener bland olika människor. Det bygger på att undersöka flera människors gener samtidigt i stället för var och en för sig, vilket sänker kostnaden. Nackdelen med detta är dock att en del av resultaten inte går att spåra till rätt individer. Därför krävs en undersökning för att avgöra hur stor del av resultaten som blir entydiga och se om det är en tillräckligt stor del i förhållande till den sänkta kostnaden. För att inte slösa på resurser innan det går att fastställa om det faktiskt är en hållbar process undersöks detta genom simuleringar i Java med Eclipse som utvecklingsmiljö.

1 Introduktion

1.1 Bakgrund

1000 Genomes Project är ett projekt som pågick mellan 2008 och 2015 där man skapade den största offentliga katalogen för mänsklig variation och genotypdata. Dessa data finns tillgängliga för allmänheten i form av VCF-filer (Variant Call Format) och används i många forskningsprojekt. Ett sätt att studera gener är med hjälp av microarrayer som fungerar som en platta med många små element som vardera representerar en position där det finns olika varianter. I dessa element skapas en reaktion som ger utslag på en av varian- terna, vilket ger information om huruvida genen finns eller inte. Metoden är vida använd idag och förhållandevis billig, men det är fortfarande kostsamt att testa stora populationer [1].

1.2 Motivering och målbeskrivning

Processen som beskrivs i bakgrunden kan göras effektivare genom att poola data från flera individer i samma microarray. Tanken med detta projekt är att undersöka hur många av individerna i microarrayen som ger entydiga och tillförlitliga svar när deras data poolas. Initialt avses grupperna innehålla åt- ta slumpmässiga individer var och i det fallet gör man sex microarrayer med fyra av dessa individer i varje. Detta sänker kostnaden med en fjärdedel. Om grupper om åtta ger en hög träffsäkerhet kan det vara intressant att använda ännu större grupper då det ytterligare sänker kostnaden.

(6)

Denna möjlighet undersöks med hjälp av tidigare nämnda data från 1000 Genomes Project. Dessa data läses in med det Java-baserade biblioteket HTSJDK som är skapat för att hantera bland annat vcf-filer [2]. De positio- nerna filtreras sedan för att matcha de som microarrayerna testar. Grupperna undersöks för att avgöra hur många individer som ger ett entydigt svar för var och en av de genetiska positioner som testas.

För varje markör avses också att undersöka dess Minor Allele Frequency (MAF) vilket är hur ofta den ovanligaste allelen förekommer i markören.

Detta är intressant då en hög MAF bör betyda att fler grupper har svårare att bli entydigt bestämda. Målet är att resultatet sedan ska kunna avläsas och exporteras till MATLAB för att utvärderas.

En ambition är att implementera API- och CLI-gränssnitt så att andra an- vändare lättare kan använda programmet. Det skulle också vara intressant att göra det möjligt att ändra parametern för storlek på grupper samt att testa specifika populationer ur den totala baserat på till exempel geografi.

2 Metod

2.1 Utförande

2.1.1 ReadVCF

Simuleringen skedde i utvecklingsmiljön Eclipse med programspråket Java.

Det Java-baserad biblioteket HTSJDK användes för att hantera VCF-filerna.

Kromosom 11 i 1000 Genomes Project fas 3 var de grunddata som användes [3]. Denna innehåller 4045628 markörer från 2504 olika individer. Dessa data filtrerades mot en textfil som innehåller 2612357 markörnamn som finns på microarrayen Infinium Omni2.5Exome-8 v1.3 BeadChip [4]. Funktionaliteten att göra detta inkapslades i en klass som heter ReadVCF, som skapades för detta projekt.

Detta gjordes genom att varje markörnamn skannades från textfilen in i en HashSet. Sedan läts en iterator gå över VCF-filen och jämföra varje VariantContext-objekts namn mot setet. HashSet valdes över

ArrayList då det gick fortare att söka i ett HashSet-objekt, testade kör- ningar med ArrayList tog över 10 timmar att slutföra. Det tog cirka 5

(7)

minuter att köra programmet med HashSet. Då den interna ordningen på namnlistan inte hade någon betydelse så fanns inget skäl att använda ArrayList.

Iteratorn var nödvändig att använda då arbetsminnet inte klarade av att konvertera hela VCF-filen i ett stycke. Vid varje överensstämmande markör- namn skrevs det VariantContext-objektet in i en ny VCF-fil som totalt fick 34749 markörer, vilka alla testade en specifik genotyp från samtliga individer.

2.1.2 RandomGroups

Slumpningen och utvärderingen gjordes i en annan klass som heter

RandomGroups. Programmet hämtade först alla individers identifikation till en lista, sedan gjordes en lista av listor som representerade grupperna.

Där slumpades identifikationerna in med hjälp av ett tidsbaserat slumpfrö medan de togs bort ur den första listan. Sedan utvärderades varje markör genom att en iterator gick igenom varje VariantContext-objekt i VCF-filen.

För varje objekt utvärderades alla grupper och resultatet adderades och skrevs in i en textfil. Markörens MAF beräknades också genom att ta fram den ovanligaste allelen, räkna hur ofta den förekom och dividera den på antal alleler. Resultatet av MAF-beräkningen skrevs in i en separat textfil. Det gjordes totalt tre simuleringar med slumpfröna 230431854001625, 129666562967335 och 130067199881048, i den ordningen.

2.1.3 Poolning

Poolningen baserades på att dela en grupp på åtta individer i två lika stora delar, tre gånger, i olika konstellationer. De åtta individerna skrevs binärt under varandra och uppdelningen i delgrupperna skedde genom att dela in nollorna och ettorna i varsinn grupp för varje kolumn. Grupperingen visas nedan i Tabell 1. I Tabell 2 visas det tydligare hur de sex delgrupperna såg ut med de åtta individerna numrerade 0-7.

(8)

Tabell 1: Denna tabell visar hur poolningen bestämdes, där varje unik färg är en delgrupp

0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1

Tabell 2: Denna tabell visar hur alla individer (numrerade 0 till 7) i vardera grupp poolades

0, 1, 2, 3 4, 5, 6, 7 0, 1, 4, 5 2, 3, 6, 7 0, 2, 4, 6 1, 3, 5, 7

Det finns totalt fyra olika kombinationer av alleler i varje markör, de kan kallas; 0|0, 1|1, 0|1 och 1|0. 1|0 och 0|1 säger i detta fall samma sak. När resultaten 1|1 eller 0|0 fås i en delgrupp betyder det att individerna i delgruppen ger entydiga resultat. För att en individ entydigt ska kunna bestämmas mås- te denne ha en homozygot genotyp och dennes genotyp måste vara likadan som de övriga genotyperna i den delgruppen (delgrupperna syns i Tabell 2).

Homozygot betyder att båda allelerna i genotypen är likadana, alltså 0|0 eller 1|1.

Om exempelvis individ nummer 7 har genotypen 1|0 och resten av individerna har genotypen 0|0 kommer det inte gå att säga något om delgrupperna med individ 7 i sig (de kommer att få resultatet 0|1 eller 1|0). Resten av delgrupperna kommer dock att vara entydiga (de kommer att få resultatet 0|0), vilket innebär att individerna 0-6 går att återskapa, men inte individ 7.

Detta utnyttjades för att analysera hur många individer som gick att enty-

(9)

digt bestämmas för respektive markör, och hur MAF påverkade resultaten.

Detta visades sedan med hjälp av plottar i MATLAB.

2.2 Verktyg

2.2.1 MATLAB

MATLAB är ett program som används för att göra matematiska beräkningar och skapa plottar utifrån resultaten. Det används i detta projekt för att utvär- dera resultaten av simuleringarna som görs i Java. Textfilerna med resultaten av utvärderingarna som görs i Java läses in i MATLAB. Sedan skapas plottar utifrån dessa data som visar hur andel individer som ger entydigt bestämda resultat förändras mellan de olika markörerna samt vid olika MAF.

2.2.2 Java

Java är ett objektorienterat programspråk, vilket innebär att program består av ett eller flera objekt som samspelar. Dessa objekt och deras egenskaper definieras i sin tur i olika klasser. En unik sak med Java är att när det språkets program kompileras skapas Javabytekod som sedan exekveras av en virituell dator kallad Java Virtual Machine (JVM). Då JVM alltid exekverar programmen gör det Javaprogram väldigt portabla då de kan köras på andra platformar utan att kompileras om [5].

2.2.3 HTSJDK

HTSJDK är ett Java-baserat bibliotek som bland annat har klasser som kan hantera VCF-filer. VCF är ett filformat som i princip är ett tabellformat i text där varje markör har data för varje individ lagrad. HTSJDK används för all inläsning och manipulering av VCF-filerna då de är för stora för att öppna på annat sätt. HTSJDK klarar också av att hantera det komprimerade vcf.gz- formatet vilket gör att ursprungsdata kunde behållas i en avsevärt mycket mindre storlek lokalt [2].

2.2.4 Eclipse

Eclipse är en utvecklingsmiljö som fungerar bra när man arbetar med pro- grammeringsspråket Java. Det är den miljön och det språket som har använts mest i detta projekt. Det finns olika vyer och perspektiv i Eclipse, vilket gör

(10)

att användaren kan anpassa miljön efter hur projektet är uppbyggt. I detta fall användes standardperspektivet Java för att skriva klasserna och bygga projektet.

2.2.5 GitHub

GitHub är en hemsida för externa Git-repositorier. Git fungerar lite annorlun- da än många andra versionhanteringsprogram i att den för varje version tar en ögonblicksbild på hela projektet istället för att bara spara ändringar i de ändrade filerna. Detta gör att man lätt kan återgå till tidigare versioner av samtliga filer som tillsammans bygger upp projektet. Det har bra integration med Eclipse och används för versionshantering och hjälper författarna att hålla varandra uppdaterade. GitHub möjligör att dela projekt med andra som då får möjlighet att föreslå ändringar i koden (push) samt uppdatera sin egen kod utifrån projektets utveckling (pull) [6].

2.2.6 Maven

Apache Maven är ett program som kan användas för att bygga mjukvaru- projekt med flera komponenter. Med M2Eclipse integrerades det direkt till Eclipse och användes för att sköta beroenden för projektet. Beroenden finns för att projektet ska hitta filer och klasser som används i projektet. Detta underlättar avsevärt arbetet då fler än en användare ska jobba med samma projekt [7].

(11)

3 Resultat

3.1 Andel entydigt bestämda individer

Andel entydigt bestämda individer i procent

0 10 20 30 40 50 60 70 80 90 100

Antal markörer

0 500 1000 1500 2000 2500 3000 3500 4000

Figur 1: Denna figur visar ett histogram över markörfördelningen i förhål- lande till andel entydigt bestämda individer för den första simuleringen, med antal markörer på y-axeln och andel entydigt bestämda individer på x-axeln.

Figur 1 visar resultatet av den första simuleringen i form av ett histogram.

Resultatet visar markörfördelningen i förhållande till andel entydigt bestäm- da individer uppdelat i intervall om 2,5% entydigt bestämda individer. Till exempel är det ungefär 4000 markörer som har mellan 2,5% och 5% entydigt bestämda individer och ungefär 1500 markörer som har mellan 10% och 12,5% entydigt bestämda individer.

(12)

1

0 10 20 30 40 50 60 70 80 90 100

Figur 2: Denna figur visar ett lådogram över hur markörerna är uppdelade i förhållande till andel entydigt bestämda individer för den första simuleringen, med andel entydigt bestämda individer på y-axeln.

Figur 2 visar motsvarande resultat som figur 1, men i form av ett lådogram i stället för ett histogram. Resultatet visar, precis som innan, markörfördel- ningen i förhållande till andel entydigt bestämda individer, men illustrerar i detta fall vilka intervall av andel entydigt bestämda individer som hör till vilka kvartiler av antal markörer.

(13)

3.2 Minor Allele Frequency (MAF)

Minor Allele Frequency (MAF) i procent

0 5 10 15 20 25 30 35 40 45 50

0 10 20 30 40 50 60 70 80 90 100

Figur 3: Denna figur visar en punktplot över hur andel entydigt bestämda individer beror av MAF för varje markör för den första simuleringen, med andel entydigt bestämda individer på y-axeln och MAF i procent på x-axeln.

Figur 3 visar ett annat resultat av samma simulering som Figur 1 och Figur 2, i form av en punktplot. Denna illustrerar hur andel entydigt bestämda individer beror av MAF för varje markör. Varje punkt motsvarar alltså en markör.

(14)

Minor Allele Frequency (MAF) i procent

0 5 10 15 20 25 30 35 40 45 50

Antal markörer

0 200 400 600 800 1000 1200 1400 1600 1800

Figur 4: Denna figur visar ett histogram över markörfördelningen i förhållan- de till MAF i den filtrerade filen, med antal markörer på y-axeln och MAF i procent på x-axeln.

Figur 4 illustrerar markörfördelningen i förhållande till MAF i den filtrerade filen. Det är en jämnare fördelning mellan staplarna i denna figur än i Figur 1, där markörfördelning i förhållande till andel entydigt bestämda individer visas, notera att MAF ligger på 23,80% i snitt.

Medelvärdena och standardavvikelserna för de olika simuleringarnas antal entydigt bestämda individer visas nedan i Tabell 3.

Tabell 3: Denna tabell visar medelvärdena och standardavvikelserna för de olika simuleringarna

Slumpfrö Medelvärde Standardavvikelse 230431854001625 934,01 757,69

129666562967335 931,66 757,51 130067199881048 934,66 756,30

Medelvärdena från Tabell 3 ligger väldigt nära varandra och motsvarar un- gefär 37% av det totala antalet individer 2504. Alltså ger ungefär 37% av

(15)

individerna ett entydigt resultat när en simulering körs med grupper om åtta individer. Standardavvikelserna från Tabell 3 ligger även de nära varandra.

4 Diskussion

4.1 Utvärdering

Det syns i Figur 1 att andel entydigt bestämda individer är ganska jämnt fördelat bland markörerna. Det är ungefär lika många markörer som har 47,5-50% entydigt bestämda individer och 90-92,5% entydigt bestämda individer. Detta gäller dock inte för låga andel entydigt bestämda individer, där antalet markörer är högre. Detta beror troligtvis på hur utvärderingen av poolerna görs, grupper kan endast få resultaten noll, fyra, fem, sex, sju samt åtta entydigt bestämda individer. Detta beror på att minst en delpool med fyra måste vara homozygota och lika för att någon individ ska bestämmas.

Markörer med hög MAF kommer således förskjuta resultaten mot noll för varje grupp.

Det syns i Tabell 3 att de olika simuleringarna har väldigt lika resultat, vilket tyder på att programmet fungerar bra, även om det egentligen krävs en bredare statistisk analys för att fastställa detta. Medelvärdet är ungefär 37% entydigt bestämda individer för samtliga simuleringar, vilket är ett resultat som anses fullgott. Detta speciellt med tanke på att MAF i genomsnitt låg på 23,80%, vilket tyder på att många av de undersökta markörerna hade tämligen vanligt förekommande varianter.

Om det är tillräckligt hög träffsäkerhet för att kunna användas i praktiken är svårt att säga, men för markörer med låg MAF bör det kunna vara ett användbart tillvägagångssätt. I Figur 3 ser vi att en stor andel av markörerna med MAF under 5% har fler än 80% entydigt bestämda individer, vilket är ett väldigt bra resultat.

Det är också värt att notera att information erhålls även för många av de individer som inte entydigt kan bestämmas. Om alla utom en individ i en viss delpool har 0|0 (vilket ibland kan bestämmas från andra delpooler) och svaret från delpoolen blir 1|0 så kan vi inte med säkerhet bestämma den sista individens alleler. Men vi vet att den minst innehåller en etta, vilket i

(16)

biologiska sammanhang kan vara viktig information i sig.

4.2 Felkällor

Ett oroväckande resultat är markörerna i nedre vänstra hörnet på Figur 3.

Vi kan inte med säkerhet förklara hur en markör får så få bestämda individer med så låg MAF. Det kan bero på att vissa markörer inte har komplett data på samtliga individer och att ett tomt anrop till dessa i utvärderingen förstör hela gruppen. Resultaten skulle också behöva statistiskt säkerställas ytterligare för att vara fullständigt pålitliga.

5 Slutsatser

5.1 Upplevelse av verktyg

5.1.1 Överblick

Då projektet startade så var båda författarna helt okunniga om några av verktygen vi arbetade med, vilket kan göra att en utvärdering av dem kan vara intressant. Java och Matlab hade vi båda erfarenhet av så dem utesluter vi ur utvärderingen då dessa fungerade som väntat. Det vi först kom i kontakt under arbetet med var Eclipse. Detta program har stora fördelar över DrJava som vi tidigare använt. Det finns mycket fler valmöjligheter och integrerade tillbehör, vilket underlättar arbetet oerhört när allt fungerade. På grund av alla valmöjligheter tar det dock avsevärt mycket längre tid att vänja sig vid än DrJava som man i princip bara startar och börjar skriva kod i. Vi använder Eclipse istället för DrJava bland annat för att DrJava inte har ett bra sätt att importera hela HTSJDK på.

5.1.2 GitHub

Det andra nya verktyget vi kom i kontakt med var Git och närmare bestämt GitHub. GitHub använder vi primärt för att dela projektet mellan förfat- tarna. Eclipse direktintegration med GitHub underlättar vid varje pull och push av projektet och är väldigt smidigt att använda när allt är konfigurerat.

Konfigureringen tar dock ganska lång tid och det är svårt att veta om man allting är rätt inställt. Vi har ingen version av vårt projekt som helt havererar

(17)

utan att vi enkelt kan spåra varför. Hade vi haft det hade tiden det tog att konfigurera git troligen känts mer väl spenderad.

5.1.3 Maven

Maven var det tredje nya verktyg vi bekantade oss med. Det var något enklare konfigurera än Git. Vårat projekt har ganska få delar så vi använde det primärt för att importera HTSJDK som ett beroende på ett väldigt enkelt sätt.

5.1.4 HTSJDK

HTSJDK var det sista obekanta verktyget vi använde. Det har exakt allt vi behöver i form av klasser för att läsa in data och bearbeta dem efter våra behov och mer därtill. Det stora hindret för att använda HTSJDK är dock inte att förstå upplägget mellan olika klasser ur ett programmeringsperspek- tiv. Det svåra för oss är att som noviser inom biologi så förstår vi inte alltid termerna de använder sig av i dokumentationen. Detta resulterar i att vi ibland måste pröva ett par olika metoder och skriva ut det de returnerar för att se om det är den metod vi söker. Att HTSJDKs dokumentation är riktad till personer med djupare kunskaper i biologi får dock anses som fullständigt rimligt.

5.2 Utveckling

Som nämnt i inledningen skulle det vara intressant att utveckla programmet så att det blir möjligt för användaren att ställa in hur många individer som ska finnas i varje grupp. Utvärderingsloopen i detta projekt blev lite för specifik och fungerar endast för grupper om åtta individer. För att lyckas med denna utveckling krävs det att det upptäcks ett samband för hur villkoren för entydighet ändras när man ökar antalet individer i varje grupp.

Det skulle också vara intressant att kunna ställa in olika parametrar, till exempel geografi, för att kunna göra simuleringar som är mer specifika för olika situationer.

Gränssnitten, även om de existerar skulle också kunna utvecklas. I nuvarande

(18)

form tar de båda programmen in argument för in- och utdata. För använda- re helt utan programmeringsbakgrund kanske ett grafiskt gränssnitt skulle kunna implementeras, detta i kombination med valmöjligheter för grupper skulle ge produkten ett professionellare intryck.

Det skulle också kunna utvecklas ett mått för den information som man får om vissa individer som inte entydigt kan bestämmas samt ett sätt att beräkna detta.

Fler körningar av programmet mot fler kromosomer skulle också behövas, samt att biologer som jobbar med microarrayer utvärderar ifall resultaten är bra nog för att metoden bör praktiseras.

6 Referenser

[1] IGSR: The International Genome Sample Resource. The 1000 Genom- es Project[Internet] Cambridgeshire: European Molecular Biology Labo- ratory; 2008[Uppdaterad 2016-05-26, citerad 2016-05-26]. Hämtad från:

http://www.1000genomes.org/about

[2] Samtools. A Java API for high-throughput sequencing data (HTS) formats.[Internet].[citerad 2016 25 maj]. Hämtad från:

https://samtools.github.io/htsjdk/

[3] IGSR: The International Genome Sample Resource. The 1000 Genomes Project[Internet] Cambridgeshire: European Molecular Biology Labora- tory; 2008[Uppdaterad 2015-05-27, citerad 2016-05-26]. Hämtad från:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr11.phase3_- shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz

[4] Illumina. Infinium Omni2.5-8 v1.3 Support Fi- les[Internet] Hämtad från: [uppdaterad 2016 mars 10;

citerad 2016 9 mars] ftp://webdata2:webdata2@ussd- ftp.illumina.com/downloads/productfiles/humanomni25/v1-

3/infiniumomni2-5-8-v1-3-a1-455-locus-report.zip

[5] Skansholm J. Java direkt med Swing. 8. uppl. Lund: Studentlitteratur AB, 2014.

(19)

[6] Chacon S, Straub B. Pro Git [Internet]. New York:

Apress; 2014 [citerad 2016 maj 25] GitHub. Hämtad från:

https://progit2.s3.amazonaws.com/en/2016-03-22-f3531/progit- en.1084.pdf

[7] The Eclipse Foundation. M2Eclipse[Internet].[uppdaterad 2015 sep 21;

citerad 2016 maj 25] Hämtad från: http://www.eclipse.org/m2e/

Appendices

A ReadVCF

/**

*@author Jonas Giek & Simon Stromstedt Hallberg

*attempt at reading a VCF

*/

package poolinggenomes;

import htsjdk.samtools.util.CloseableIterator;

import htsjdk.variant.vcf.VCFFileReader;

import htsjdk.variant.variantcontext.VariantContext;

import htsjdk.variant.variantcontext.writer.VariantContextWriter;

import htsjdk.variant.variantcontext.writer.VariantContextWriterBuilder;

import htsjdk.variant.variantcontext.writer.Options;

import htsjdk.variant.vcf.VCFHeader;

import htsjdk.samtools.SAMSequenceDictionary;

import java.util.HashSet;

import java.io.*;

import java.util.Scanner;

public class ReadVCF { /**

* This is a program to filter an existing .vcf-file against a textfile

* containing the names of the variants you want to keep.

* @param args The first argument is the filename of marker names that the program will filter the .vcf-file against.

* @param args The second argument is the filename of the .vcf or vcf.gz-file you want to filter.

* @param args The third argument is the name of the output file .

*/

public static void main(String[] args) throws

(20)

FileNotFoundException {

Scanner scan = new Scanner(new File(args[0]));

HashSet<String> names = new HashSet<String>();

int countfilter= 0;

int countgrundfil= 0;

while (scan.hasNext()){

scan.nextInt();

names.add(scan.next());

countfilter++;

scan.nextLine();

}

scan.close();

File vcfFile = new File(args[1]);

VCFFileReader reader = new VCFFileReader(vcfFile);

SAMSequenceDictionary dic = VCFFileReader.

getSequenceDictionary(vcfFile);

VCFHeader header = reader.getFileHeader();

VariantContextWriterBuilder builder = new VariantContextWriterBuilder()

.setReferenceDictionary(dic)

.setOption(Options.INDEX_ON_THE_FLY);

VariantContextWriter writer = builder .setOutputFile(args[3]) .build();

CloseableIterator<VariantContext> iter = reader.iterator ();

VariantContext variant = null;

writer.writeHeader(header);

int count= 0;

while (iter.hasNext()) { variant = iter.next();

if (names.contains(variant.getID()) && variant!=

null) {

writer.add(variant);

count++;

}

countgrundfil++;

}

System.out.println(count);

System.out.println(countfilter);

System.out.println(countgrundfil);

System.out.println();

System.out.println(variant.getID());

reader.close();

} }

(21)

B RandomGroups

/**

* @author Jonas Giek & Simon Stromstedt Hallberg

*///

package poolinggenomes;

import htsjdk.samtools.util.CloseableIterator;

import htsjdk.variant.vcf.VCFFileReader;

import htsjdk.variant.variantcontext.VariantContext;

import htsjdk.variant.variantcontext.Allele;

import java.io.*;

import java.util.Set;

import java.util.ArrayList;

import java.util.List;

import java.util.Random;

class RandomGroups { /**

* This a program to calculate how many of the indivduals in a certain vcf-file

* could get their genotype correctly determined with randomised pooling in microarrays.

* Output files contain how many correctly determined individuals each marker got as well

* as the Minor Allele Frequency of the marker.

* @param args The argument is the input file that will be evaluated.

*/

public static void main(String[] args) throws FileNotFoundException {

long seed = System.nanoTime();

// long seed = null;

File vcfFile = new File(args[0]);

VCFFileReader reader = new VCFFileReader(vcfFile, false);

CloseableIterator<VariantContext> iter = reader.iterator();

VariantContext variantForNames = iter.next();

Set<String> names = variantForNames.getSampleNames();

ArrayList <String> namelist = new ArrayList(names);

int size = namelist.size();

int shrinkingSize = namelist.size();

int sizeOfGroups = 8;

Random rdm = new Random(seed);

String seedToNameCount =

"markcount"+String.valueOf(seed)+".txt";

String seedToNameMAF = "maf"+String.valueOf(seed)+".txt";

ArrayList <ArrayList<String>> groups =

(22)

new ArrayList<ArrayList<String>>();

for (int i = 0; i < (size/sizeOfGroups); i++) {

ArrayList<String> L = new ArrayList<String>();

for (int j = 0; j <sizeOfGroups; j++) { int ind = rdm.nextInt(shrinkingSize);

String S = namelist.get(ind).toString();

L.add(S);

namelist.remove(ind);

shrinkingSize--;

}

groups.add(L);

}

iter.close();

CloseableIterator<VariantContext> iter2 = reader.iterator();

VariantContext variant = null;

int countmarkers=0;

double maf;

List<Allele> alleles = null;

PrintWriter outputStream=new PrintWriter(seedToNameCount);

PrintWriter outputStream1=new PrintWriter(seedToNameMAF);

//PrintWriter outputStream=new

PrintWriter(seedToName+"replica.txt");

int count;

while (iter2.hasNext()) { variant = iter2.next();

alleles = variant.getAlleles();

variant.getCalledChrCount(alleles.get(0));

if (variant.getCalledChrCount(alleles.get(0))

>=variant.getCalledChrCount(alleles.get(1))) { maf = (double) variant.getCalledChrCount (alleles.get(1))/variant.getCalledChrCount();

} else {

maf = (double) variant.getCalledChrCount (alleles.get(0))/variant.getCalledChrCount();

}

outputStream1.println(maf);

countmarkers++;

count=0;

for (int j = 0; j < groups.size(); j++) {

if (variant.getGenotype(groups.get(j).get(0)).isHom()

&& (variant.getGenotype(groups.get(j).get(0)).

sameGenotype(variant.getGenotype(groups.get(j).get(1)))

&& variant.getGenotype(groups.get(j).get(0)).

|| variant.getGenotype(groups.get(j).get(0)).

(23)

sameGenotype(variant.getGenotype(groups.get(j).get(6))))) count++;

(24)

&& variant.getGenotype(groups.get(j).

get(3)).sameGenotype(variant.getGenotype(groups.get(j).get(2)))

|| variant.getGenotype(groups.get(j).

get(3)).sameGenotype(variant.getGenotype(groups.get(j).get(7))))) count++;

&& (variant.getGenotype(groups.get(j).

(25)

(26)

sameGenotype(variant.getGenotype(groups.get(j).get(3)))&&

variant.getGenotype(groups.get(j).get(7)).

}

outputStream.println(count);

}

outputStream.close();

outputStream1.close();

System.out.println(seed);

System.out.println(countmarkers);

reader.close();

} }