Page ThesisCNVs

Detektion von Copy Number Variationen in Sequenzierdaten

Student

Kerstin Neubert

Betreuer

Anne-Katrin Emde, Prof. Dr. Knut Reinert

Zusammenfassung

In der Masterarbeit soll eine Methode zur Detektion von Copy number variations (CNVs) in Second Generation Sequencing Daten implementiert werden. Der entwickelte Algorithmus wird anschließend auf simulierte und reale Sequenzierdaten von einem Tumorgenom und Daten aus dem 1000 Genomes Projekt angewendet.

Beschreibung

Second-generation sequencing Technologien ermöglichen es, die Positionen von single nucleotide polymorphisms (SNPs) und strukturellen Variationen wie z.B. copy number variations (CNVs) im Genom mit einer größeren Genauigkeit im Vergleich zu arraybasierten Methoden wie z.B. SNP-Microarrays und der whole-genome array comparative genome hybridization (aCGH) zu bestimmen. Jedoch erfordert die Analyse der dabei anfallenden großen Datenmengen neue effiziente Algorithmen. In dieser Arbeit soll ein Algorithmus für die Analyse von CNVs in second-generation sequencing Daten entwickelt und in C++ implementiert werden. Dabei werden folgende Teilaufgaben bearbeitet:

1. Entwicklung und Implementierung eines Algorithmus zur Detektion von CNVs in 2nd-generation sequencing Daten. Fragestellungen tun sich auf bezüglich der Wahl der Fensterbreite (dynamisch/statisch) und hinsichtlich der unterschiedlich guten "Sequenzierbarkeit" sowie "Mappbarkeit" von bestimmten Sequenzen.
2. Bestimmung der Spezifität und Sensitivität des entwickelten Algorithmus anhand simulierter Daten, bei denen die Position von CNVs bekannt ist
3. Anwendung des Algorithmus zur Detektion von CNVs in verschiedenen 2nd-generation sequencing Daten (ein Tumorgenom und Daten aus dem 1000 Genomes Projekt)
4. Bewertung des Algorithmus (Laufzeit, Vergleich mit anderen Algorithmen)

Literatur

[1] Paul Medvedev, Monica Stanciu, and Michael Brudno: Computational methods for discovering structural variation with next-generation sequencing. Nature Methods 2009

[2] Seungtai Yoon, Zhenyu Xuan, Vladimir Makarov, Kenny Ye, and Jonathan Sebat: Sensitive and accurate detection of copy number variants using read depth of coverage. Genome Research 2009

[3] Korbel et. al.: Paired-end Mapping Reveals Extensive Structural Variation in the Human Genome. Science 2007

Arbeitsplan

10.12.-18.12. Organisatorisches, Literaturrecherche

19.12.-21.02. Entwicklung des Algorithmus und Implementierung in C++, Testläufe mit simulierten Daten

22.02.-14.03. Anwendung auf reale Daten (z.B. 1000Genomes Projekt), Optimierung bzgl. Laufzeit

15.03.-30.04. Bewertung der Ergebnisse (Sensitivität, Spezifität, Laufzeit, Vergleich mit anderen Ansätzen, Datenbankeinträgen für CNVs)

1.05.-23.05. Masterarbeit schreiben

24.05.-9.06. Korrekturen

10.06. Abgabe

Wöchentliche Berichte

[4.01.-10.01.]

• Suche nach verfuegbaren CNV-Tools fuer single-end Daten
• Literatur

[11.01.-17.01]

• Testlauf von CNV-seq: Detektion von copy number Veraenderungen in zwei NGS-Datensaetzen (SOLiD)
• 15179 CNVs im humanen Genom, 1,5% (hg19), (Mittelwert 3 kb, max. 231 kb)

[18.01.-24.01.] / [25.01- 31.01.]

Vergleich verschiedener Ansaetze zur CNV-Detektion basierend auf read depth:

• CNV-seq [Xie and Tammi, Bioinformatics 2009]
• Event-wise testing [Yoon et. al., Genome Res. 2009]
• SegSeq [Chiang et. al., Nat. Methods 2009]
• Applied Biosystems CNV Tool [http://solidsoftwaretools.com/gf/project/cnv/]

[1.-7.02.]

• Seqan Einarbeitung
• statistische Verteilungen und Tests mit boost C++ libraries

[8.-14.02.]

Algorithmus [CNV-Detektion in zwei NGS-Datensaetzen]

Input: mapped reads, parameter (log2_threshold, p_threshold)

Output: CNVs (Start-Ende, copy number)

1. Zaehle reads in Fenstern entlang des Genoms in beiden Proben (Startpositionen): n_x, n_y

• feste Festergroesse, nicht-ueberlappend
• feste Fenstergroesse, ueberlappend (sliding windows)
• dynamische Fenstergroesse, "local change points" (Teststatistik) oder konstante Anzahl "uniquely mapped" Positionen im Referenz-NGS-Datensatz

2. Berechne das Verhaeltnis der Anzahl an reads in den Fenstern der beiden Datensaetze: r=log2(n_x/n_y)

3. Auflistung von Festern mit r >= |log2_treshold| und p <= p_threshold

4. iterative Zusammenfassung benachbarter Fenster zu CNV-Regionen, falls |r| gleich ist

Algorithmus [CNV-Detektion in einem NGS-Datensatz]

wie bei 2 Datensaetzen, ein Datensatz wird dabei simuliert

Simulation von single-end Daten:

• stat. Verteilung entsprechend Original-Sequenzierdaten (Poisson?)
• coverage
• read size
• sequencing bias (GC)

[15.02.-21.02.]

• Einlesen von Read-Alignments (SAM-Format) in FragmentStore
• BAM-Format: Umwandlung von BAM- in SAM-Format mit SAMTools
• GFF-Format: SOLiD gff2sam Skript

[22.02.-28.02.]

• Programmablauf:

     read (CharString file, FragmentStore &me) {
        // parse aligned reads into fragment store from SAM format
     }
  
     initWindows (FragmentStore &me, TWindowStore &windowStore, bool overlap) {
        // identify (non) overlapping windows of dynamic/constant size
     }
  
     countReads (FragmentStore &me, TWindowStore &windowStore, bool overlap) {
        // count read start positions in windows for each contigID
     }
  
     mergeWindows (TWindowStore &windowStore, TWindowStore &merged) {
        // Identify overlapping/adjacent windows with equal read count ratio
     }
  
     callCNVs (TCNVStore &cnvs, TWindowStore &merged) {
        // call CNVs for overlapping/adjacent windows (position, copy number)
     }
  
     filterCNVs (TCNVStore &cnvs, unsigned size, double pval) {
        // Filter CNVs by user defined filtering criteria (size, pvalue)
     }
  
     writeCNVs (CharString file, cnvStore &cnvs) {
        // write cnv results in tab separated file
     }
     

[1.03.-7.03.]

• Definition von Klassen für Fenster (konstante/dynmische Größe)
• WindowStore für alle Fenster
• Testdatensatz von je 1 SOLiD Run von Tumor/Kontrolle für Chromosom 9 (mapping auf hg18) erzeugt (ca. 2 Mill. Reads)

[8.03.-14.03.]

• Implementierung der countReads-Funtion zum Zählen der Reads in Fenstern entlang des Genoms

[15.03.-21.03.]

• Klasse CNVoptions für Parameter, Funktionsmoden
• Normalisierung der Read counts mit GC nach countReads im Programmablauf eingefuegt

[22.03.-28.03.]

• Implementierung der GC-Normalisierung der read counts

[29.03.-04.04.]

• Berechnung eines lokalen CNV Scores fuer einzelne Fenster (Z-score)
• z = ( rc-mean(rc) ) / sd(rc)

[05.04.-11.04.]

• 1.Testdatensatz (SOLiD) von Chromosom 9 auf 15 Test- und 18 Kontrol- SOLidruns (coverage ca. 20x) erweitert
• 2. Testdatensatz: Chromosom 1 vom CEU-Trio (1000 Genomes-Projekt, Solexa) mit SAMtools von BAM in SAM-format umgewandelt --> Fehler beim Einlesen der SAM-Datei in FragmentStore
• Umwandlung des Z-scores der Fenster in upper-tail probability P(Z >= z_i) für CNVs mit erhöhter Kopiezahl (Gains)
• Implementierung von Event-wise testing (EWT) für CNV events in aufeinanderfolgenden Fenster

[12.04.-18.04.]

• Einlesefunktion für SAM-Dateien überladen (für Contig-ID, Anfangsposition, Endposition)
• Einlesefunktion für FASTA-Sequenzen in ContigStore überladen, so dass nur die Sequenz eines Contigs eingelesen wird
• Speicherbelastung bei realen Daten (mit coverage >20x) sehr groß → stückweises Einlesen in FragmentStore für sequentielle Analyse (z.B. Reads in Fenstern von 1 Million Basen auf einmal einlesen und dann auswerten)

[19.04.-25.04.]

• Funktionen zum Erzeugen von log-Dateien/Tabellen für Testläufe: Statistik des Input-Datensatzes (coverage, GC% im Genom), Read counts (Rohdaten und normalisiert), Scores, p-Werte, Matrix für CNV events
• mapping quality der Reads einbezogen (minimale mapping quality, Berechnung der duchschnittlichen mapping quality der Reads in Fenstern)
• Funktion zum stückweisen Einlesen der Reads (mit maximal auf einmal einzulesende Reads als Parameter, z.B. 1 Million) geändert
• Fehler beim Einlesen der SAM-Dateien des 1000Genomes-Datensatzes beseitigt durch auskommentieren von "insert alignmentgaps" in eigener Einlesefunktion

[26.04.-2.05.]

• Zusammenlegen sich überschneidender oder benachbarter Copy number events für jeweils Deletion und Duplikation
• Ausgabe der 'mergedEvents' mit duchschnittlichem Read count und p-Wert in Datei

[3.05.-09.05.]

• Mappability-track für 50 bp single-end mit Bowtie (wiggle-format) auf hg18 erstellt
• Simulation von CNVs auf zufälliger Sequenz (Einfügen von Duplikationen in Sequenz, Simulation von Reads mit ReadSimulator, Mapping mit Bowtie)

[10.05.-16.05.]

• variable Fenstergröße basierend auf Mappability
• Normalisierung der read counts mit Mappability

[17.05-23.05.]

• t-test (für 2 unabhängige Stichproben) auf Mittelwerten der read counts für Filterung von "dimorphic" Copy number events in 2 Datensätzen implementiert
• Filterung der CNVs basierend auf der relativen Abweichung der median read counts des CN events vom Mittelwert der read counts im gesamten Genom und dem p-Wert des CN events (Filterparameter)
• Ausgabe der Ergebnisse für CNVs in gff-Datei

[24.05.-30.05.]

• Teslauf auf zwei 1000Genomes-Datensätzen (Chromosom 1, 180/156 Millionen reads) mit 100bp-Fenstern → Laufzeit: ca. 10 h
• Einlesen der Daten in FragmentStore ist der Flaschenhals bzgl. Laufzeit des gesamten Programms
• Berechnung der "mappability" für das Genom abhängig davon, wie "uniqueness" einer genom. Position definiert wird, wie der mapping Algorithmus mit Mismatches (Edit/Hamming-Distanz) und mit Multireads (Reads, die an mehreren Stellen im Genom passen) umgeht
• falls eine Normalisierung bzgl. "uniquely mappable positions" angewendet werden soll, so muss der Benutzer des CNV-Tools selbst den Mappability-Track mit denselben Parametern (Anzahl Mismatches) erstellen, die beim Mapping der Sequenzdaten verwendet wurden

[31.05.-06.06.]

• Implementierung eines CNV-Simulators für die künstliche Erzeugung von CNVs (Duplikationen/Deletionen mit geg. Länge) auf einer geg. Sequenz mit Parametern: Position, Länge, CN (copy number), Fehlerstatistik für Sequenzierfehler

[07.06.-13.06.]

• Erweiterung des CNV-Simulator für die Simulation von beliebig vielen CNVs eines Typs (gleiche Länge, Kopiezahl, zufällige Position) in einer gegebenen Sequenz und Simulation von single-end reads mit Sequenzierfehlern (nur Mismatches) fertiggestellt

[14.06.-20.06.]

• Simulation von 1000 Duplikationen (Länge 5 kb) mit zufälligen Positionen auf Chromosom 1, Simulation und Mapping (mit Bowtie - Parameter: -f -y -k 2 -m 1 -v 3 --best --sam) von single-end reads auf manipulierter und Originalsequenz (max. 3 Fehler), coverage ca. 2x
• Ergebnis mit simulierten Daten: 90% TP, 10% FP CNV events, ausserdem 101 CNVs nicht detektiert
• Problem bei Simulation auf echter Sequenz: Gaps, z.B. Centromer (sollten ausgelassen werden)

[21.06.-27.06.]

• Simulation von 200 Duplikationen (Länge 1 kb) mit zufälligen Positionen auf zufälliger Sequenz (10 Mill. bp), Simulation und Mapping von SE-reads auf manipulierter und ursprünglicher Sequenz (max. 3 Fehler) mit Razers (Parameter -i 94 --unique -of 4) → coverage ca. 8x
• Ergebnis: 203 events, davon 197 TP, 6 FP, ausserdem 3 CNVs nicht gefunden (3 FN)

[28.06.-04.07.]

• CNV-Simulator: Simulation von CNVs auf echter Sequenz (Chromosom 20), nicht in repeat-maskierten Bereichen (Centromer etc.), mit Sequenzierfehler 0.01, gemischte Simulation (Duplikationen, Deletionen mit beliebiger Grösse) möglich

[05.07.-11.07.]

• Chi-Quadrat Test für Test, ob Standardabweichung der read counts in Test- und Kontrolldatensatz gleich sind (je nachdem t-test für Stichproben mit gleicher oder unterschiedlicher Standardabweichung)
• Neu-Implementierung der appendReads()-Funktion für Einlesen der Reads in FragmentStore funktioniert nicht ohne Sequenzen (werden eigentlich nicht benötigt)

[12.07.-18.07.]

• Korrektur von Filterfunktion: 2-sided t-test durch one-sided t-test ersetzt (bei 2-sided zu strenger Filter)
• Erweiterung des CNV-Programms mit weiteren Algorithmen ermöglicht, indem die Funktion für die Detektion der CNVs mit Tag aufgerufen wird: callCNVs( WindowStore, EventStore, Options, SampleId, Mode ), in EventStore sind die Ergebnisse gespeichert, im WindowStore die read counts der einzelnen Fenster

[19.07.-25.07.]

• Vergleichsprogramme für CNV-Detektion getestet:

1. cnv-seq (Fenster-basierter Ansatz basierend auf log-Ratios der read counts von Test- und Kontrolldatensatz)

2. DNAcopy (Circular binary segmentation; da ursprünglich auf Microarrays angewendet, müssen die logRatio-Inputdaten erst aus SAM-Alignment-Dateien berechnet werden)

[26.07.-01.08.]

• weiteres Vergleichsprogramme für CNV-Detektion getestet:

3. SOLiD-CNV-TOOL (implementiert ein HMM, funktioniert nur auf hg18, weil binäre mappability-Dateien nur für hg18 im Programmpaket enthalten)