fn main() {
    let none_value: Option<i32> = None;

    println!("{}", none_value.expect("Option war None"));
}

Sicherere Alternativen zu unwrap

match.

Siehe hierzu Kapitel01-06.

fn main() {
    let some_value = Some(42);

    match some_value {
        Some(value) => println!("Wert: {}", value),
        None => println!("Kein Wert gefunden"),
    }
}

unwrap_or

Falls ein Wert nicht existiert, liefert unwrap_or einen Standardwert zurück.

fn main() {
    let none_value: Option<i32> = None;
    println!("{}", none_value.unwrap_or(0)); // Gibt 0 aus.
}

unwrap_or_else

Führt eine Funktion aus, um einen Standardwert zu berechnen, wenn der Wert nicht existiert.

fn main() {
    let none_value: Option<i32> = None;
    println!("{}", none_value.unwrap_or_else(|| 42)); // Gibt 42 aus.
}

?-Operator

Der ?-Operator kann in Funktionen verwendet werden, die einen Result zurückgeben. Er kann verwendet werden, um den Wert aus einem Result oder Option zu extrahieren, ohne jedes Mal explizit match schreiben zu müssen. Wenn ein Fehler auftritt, gibt der ?-Operator den Fehler direkt aus der umgebenden Funktion zurück.

fn divide(a: f64, b: f64) -> Result<f64, String> {
    if b == 0.0 {
        Err(String::from("Division durch Null"))
    } else {
        Ok(a / b)
    }
}

fn calculate() -> Result<(), String> {
    let result = divide(10.0, 2.0)?; // Extrahiert den Wert oder gibt den Fehler (gleich) zurück.
    println!("Ergebnis: {}", result);
    Ok(()) // Rückgabe im Erfolgsfall
}

fn main() {
    match calculate() {
        Ok(()) => println!("Berechnung erfolgreich!"),
        Err(e) => println!("Fehler: {}", e),
    }
}

Interfaces (Schnittstellen)...?

In Sprachen wie Java definiert eine Interface (Schnittstelle) die Signaturen von Methoden (ohne Implementierung), die von einer implementierenden Klasse bereitgestellt werden müssen.

public interface Animal {
    void makeSound();
}

public class Dog implements Animal {
    @Override
    public void makeSound() {
        System.out.println("Wuff!");
    }
}

public class Cat implements Animal {
    @Override
    public void makeSound() {
        System.out.println("Miau!");
    }
}

public class Main {
    public static void main(String[] args) {
        Animal dog = new Dog();
        Animal cat = new Cat();

        dog.makeSound(); // Ausgabe: Wuff!
        cat.makeSound(); // Ausgabe: Miau!
    }
}

In Python gibt es keine Schnittstellen im engeren Sinne, jedoch kann man mit verschiedenen Konstrukten Ähntliches implementieren (vgl Duck typing).

Traits

In Rust haben Traits in etwa den Stellenwert von Interfaces, jedoch werden sie anders implementiert.

Man kann bestehende Traits für eigene Datentypen implementieren oder aber neue Traits definieren. Hier wollen wir den Trait FromString aus der Standardbibliothek für eine eigene Struktur implementieren.

Der Trait FromString ist von der Standardbibliothek.

#![allow(unused)]
fn main() {
pub trait FromStr: Sized {
    type Err;
    fn from_str(s: &str) -> Result<Self, Self::Err>;
}
}

Die parse-Method erwartet nur, dass ein Datentyp diesen Trait implementiert.

#![allow(unused)]
fn main() {
pub fn parse<F: FromStr>(&self) -> Result<F, F::Err> {
    FromStr::from_str(self)
}
}

Stellen wir uns folgende Struktur vor

#![allow(unused)]
fn main() {
struct DnaString {
    dna: String;
}
}

Um unsere Struktur mit parse verwenden zu können, müssen wir FromStr implementieren. Unser Code überprüft, ob alle Zeichen in dem zu parsenden String A, C, G, order T sind.

#![allow(unused)]
fn main() {
use std::str::FromStr;

struct DnaString {
    pub dna: String,
}

#[derive(Debug)]
struct DnaParseError;

impl FromStr for DnaString {
    type Err = DnaParseError;

    fn from_str(s: &str) -> Result<DnaString, DnaParseError> {
        if s.chars().all(|c| matches!(c, 'A' | 'C' | 'G' | 'T')) {
            Ok(DnaString { dna: s.to_string()})
        } else {
            Err(DnaParseError)
        } 
    }
}

#[cfg(test)]
mod test {
    use super::DnaString;

    #[test]
    fn test_dna_parse() {
        let valid_dna = "ATGCTCGCTAG";
        let dna_struct = valid_dna.parse::<DnaString>();
        assert!(dna_struct.is_ok());
        let dna_struct = dna_struct.expect("Could not unwrap DNA");
        assert_eq!(valid_dna, dna_struct.dna);
    }

    #[test]
    fn test_invalid_dna_parse() {
        let valid_dna = "ATGCT??CGCTAG";
        let dna_struct = valid_dna.parse::<DnaString>();
        assert!(dna_struct.is_err());
    }
}
}

Übung 1

Der folgende Code kompiliert nicht

#![allow(unused)]
fn main() {
let d = DnaString { dna: "TGCAACGT".to_string()};
writeln!("{}", d);
}

Verwenden Sie die Fehlermeldung des Rustkompilierers um herauszufinden, warum der Code nicht funktioniert. Es gibt mindests zwei Korrekturen - finden Sie sie und schreiben Sie entsprechenden Rustcode.

Sortieralgorithmen

Wir werden als Übung mehrere Sortieralgorithmen in Rust implementieren. Sinn und Zweck der Überung ist nicht unbedingt, Sortieralgorithmen zu erlernen, sondern zu lernen, wie man mit Traits in Rust umgeht und wie der Debugger eingesetzt werden kann.

Bevor wir mit diesem Abschnitt beginnen, muss der Debugger (CodeLLDB) als VS Code extension installiert sein, s. dazu Extensions.

Bubblesort

Bubblesort ist ein einfacher Sortieralgorithmus, der Elemente in einer Liste durch wiederholtes Vergleichen benachbarter Elemente sortiert. Dabei werden zwei benachbarte Elemente vertauscht, wenn sie in der falschen Reihenfolge stehen. Der Algorithmus wird so lange wiederholt, bis die Liste vollständig sortiert ist. Der Name "Bubblesort" kommt daher, dass kleinere Elemente in der Liste "nach oben steigen", ähnlich wie Luftblasen (bubbles) im Wasser.

Algorithmus

1. Starte am Anfang der Liste.
2. Vergleiche Element i (e_i) mit Element i+1 (e_i+1)
3. Wenn e_i größer als e_i+1 ist, vertausche sie. Andernfalls lasse sie unverändert.
4. Gehe so durch die gesamte Liste.
5. Wiederhole den Vorgang für die gesamte Liste, bis keine Vertauschungen mehr notwendig sind (die Liste ist dann sortiert).

Dieser Algorithmus hat eine Zeitkomplexität von O(n²).

#![allow(unused)]
fn main() {
pub fn bubble_sort<T: PartialOrd>(v: &mut [T]) {
    for i in 0..v.len() {
        let mut sorted = true;
        for j in 0..(v.len()-1) - i {
            if v[i] > v[i+1] {
                v.swap(i, i+1);
                sorted = false;
            }
        }
        if sorted {
            return;
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_bubble_sort() {
        let mut v = vec![4, 6, 77, 1, 13, 24];
        bubble_sort(&mut v);
        assert_eq!(v, vec![1, 4, 6, 13, 24, 77]);
    }
}
}

Übungen

Erzeuge eine Rust-Applikation namens ch2 mit cargo.

1. Führe den oben stehenden Rustcode aus.
2. Definiere einen struct-Datentype mit zwei Ganzzahl-Feldern. Konsulitere die Dokumentation für den PartialOrd-Trait (https://doc.rust-lang.org/std/cmp/trait.PartialOrd.html). Implementiere den Trait für die Struktur, so dass structs gemäß der Summe beider Felder sortiert werden. Sortiere drei Beispiel structs.
3. Schreibe einen Test der bubble-Sort-Methode für Vektoren von Ganzzahlen
4. Schreibe einen Test für die oben definierte Struktur
5. Schreibe eine Testmethode, welche mehrere Tests auf einmal durchführt.

Quicksort

• Texte teilweise von Wikipedia angepasst.
• Einführung auf YouTube.

Quicksort ist ein schneller, rekursiver Sortieralgorithmus, der nach dem Prinzip Teile und herrsche arbeitet. Im Durchschnitt führt der Quicksort-Algorithmus \(O(n\cdot log(n))\) Vergleiche durch.

Zunächst wird die zu sortierende Liste in zwei Teillisten („linke“ und „rechte“ Teilliste) getrennt. Dazu wählt Quicksort ein sogenanntes Pivotelement aus der Liste aus. Alle Elemente, die kleiner als das Pivotelement sind, kommen in die linke Teilliste, und alle, die größer sind, in die rechte Teilliste. Die Elemente, die gleich dem Pivotelement sind, können sich beliebig auf die Teillisten verteilen. Nach der Aufteilung sind die Elemente der linken Liste kleiner oder gleich den Elementen der rechten Liste.

Anschließend muss man also noch jede Teilliste in sich sortieren, um die Sortierung zu vollenden. Dazu wird der Quicksort-Algorithmus jeweils auf der linken und auf der rechten Teilliste ausgeführt. Jede Teilliste wird dann wieder in zwei Teillisten aufgeteilt und auf diese jeweils wieder der Quicksort-Algorithmus angewandt, und so weiter. Diese Selbstaufrufe werden als Rekursion bezeichnet. Wenn eine Teilliste der Länge eins oder null auftritt, so ist diese bereits sortiert und es erfolgt der Abbruch der Rekursion.

Wir haben hier den Java-Code vom oben genannten YouTube-Video leicht angepasst bzw. kommentiert wiedergegeben. Sie können optional das Programm laufen lassen, indem Sie den Code in eine Datei namens Quick.java kopieren und diese Datei mit $ javac Quick.javakompliieren und das Ergebnis mit $ java Quickausführen.

import java.util.Random;

public class Quicksort {

  public static void main(String[] args) {
    int[] numbers = [4, 6, 1, 8, 13, 33, 2];
    printArray(numbers);
    quicksort(numbers);
    printArray(numbers);
  }

  private static void quicksort(int[] array, int lowIndex, int highIndex) {
    // 1. Basisfall der Rekursion: jedes Feld mit nur einem Element ist schon sortiert!
    if (lowIndex >= highIndex) {
      return;
    }
    // 2. Das letzte ElementWir dient als Pivotelement
    int pivot = array[highIndex];
    // 3. Die zu sortierende Liste in zwei Teillisten wird in zwei Unterlisten getrennt.
    // Alle Elemente, die kleiner als das Pivotelement sind, kommen in die linke Teilliste, 
    // und alle, die größer sind, in die rechte Teilliste.
    // Nach Ausruf von 'partition' ist die Position des Pivotelements korrekt, aber die
    // linke und rechte Unterlisten sind nicht unbedingt richtig sortiert
    int leftPointer = partition(array, lowIndex, highIndex, pivot);
    // 4. Sortierte die linke Unterliste (rekursiv)
    quicksort(array, lowIndex, leftPointer - 1);
    // 5. Sortierte die rechte Unterliste (rekursiv)
    quicksort(array, leftPointer + 1, highIndex);

  }

  private static int partition(int[] array, int lowIndex, int highIndex, int pivot) {
    int leftPointer = lowIndex;
    int rightPointer = highIndex - 1;

    while (leftPointer < rightPointer) {
        // 1. Von links gehen, bis eine Zahl gefunden wird, die größer als der Pivot ist, 
        // oder der linke Zeiger mit dem rechten identisch ist
        while (array[leftPointer] <= pivot && leftPointer < rightPointer) {
            leftPointer++;
        }
        // 2. Von rechts gehen, bis eine Zahl gefunden wird, die kleiner als der Pivot ist, 
        // oder der rechte Zeiger mit dem linken identisch ist.
        while (array[rightPointer] >= pivot && leftPointer < rightPointer) {
        rightPointer--;
        }
        //  3. Die Werte vertauschen, auf die der linke und der rechte Zeiger zeigen
        swap(array, leftPointer, rightPointer);
    }
    // 4. An dieser Stelle muss das Pivotelement mit dem Wert am linken Zeiger vertauscht werden
    // (highIndex ist der Index des Pivotelements)
    if(array[leftPointer] > array[highIndex]) {
      swap(array, leftPointer, highIndex);
    }
    else {
      leftPointer = highIndex;
    }
    
    return leftPointer;
  }

  private static void swap(int[] array, int index1, int index2) {
    int temp = array[index1];
    array[index1] = array[index2];
    array[index2] = temp;
  }

  private static void printArray(int[] numbers) {
    for (int i = 0; i < numbers.length; i++) {
      System.out.println(numbers[i]);
    }
  }
}

FASTQ

FASTQ ist ein Dateiformat für den Austausch von Sequenzierungsdaten, das sowohl die Sequenz als auch eine entsprechende Bewertung der Basenqualität enthält.¹ Sein Name leitet sich vom ehrwürdigen FASTA-Format ab, einem textbasierten Format zur Darstellung von DNA- und Aminosäuresequenzen, das ursprünglich 1985 zusammen mit einem Softwarepaket für das Alignment von Nukleinsäuren und Proteinen veröffentlicht wurde Lipman et al. 1985. Das FASTA-Format ist eines der in der Bioinformatik am weitesten verbreiteten Formate überhaupt. Die erste Zeile beginnt mit einem > (Größer-als)-Symbol, direkt gefolgt von einem Namen oder einem eindeutigen Sequenzbezeichner, optional gefolgt von einer Beschreibung der Sequenz. Die folgenden Zeilen enthalten die Sequenz (in der Regel mit 70 oder 80 Zeichen pro Zeile). Eine Datei kann mehrere Sequenzeinträge enthalten, aber jeder muss mit einer >-Zeile beginnen. Die FASTA-Datei für die Homo sapiens Fibrillin 1 (\textit{FBN1}) messenger RNA beginnt zum Beispiel so

>NM_000138.4 Homo sapiens fibrillin 1 (FBN1), mRNA
AGTATTTCTCTCGCGAGAAACCGCTGCGCGGACGATACTTGAAGAGGTGGGGAAAGGAGGGGGCTGCGGG
AGCCGCGGCAGAGACTGTGGGTGCCACAAGCGGACAGGAGCCACAGCTGGGACAGCTGCGAGCGGAGCCG
AGCAGTGGCTGTAGCGGCCACGACTGGGAGCAGCCGCCGCCGCCTCCTCGGGAGTCGGAGCCGCCGCTTC
(weitere 165 Zeilen)

Das FASTQ-Format ist eine Erweiterung des FASTA-Formats, in dem zusätzlich ein numerischer Qualitätswert für jedes Nukleotid in einer Sequenz gespeichert wird. Jeder Read in der Datei wird durch vier Zeilen mit den folgenden Informationen dargestellt (ohne Längenbegrenzung):

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 

Die @-Zeile umfasst ein freies Formatfeld, das in erster Linie für den Titel oder die Kennung der gemeldeten Sequenz verwendet wird. Illumina hat eine typische Namenskonvention für die Titelzeile, die wir weiter unten erläutern werden. Die zweite Zeile enthält die aufgerufenen Basen, d. h. die Sequenz des Read. Es folgt die +-Zeile Obwohl frühe Versionen des FASTQ-Formats in dieser Zeile eine Wiederholung der Titelzeile enthielten, ist es inzwischen Standard, dass die Zeile nur das Pluszeichen (+) enthält, was die Dateigröße erheblich reduziert. Die vierte Zeile schließlich enthält die ASCII-kodierten Phred-Qualitätsbewertungen und muss die gleiche Länge wie die Sequenzzeile haben. Somit wird für jede der in der zweiten Zeile angegebenen Basen eine Qualitätsbewertung in der entsprechenden Position der vierten Zeile angegeben. In einigen früheren Versionen von FASTQ wurden verschiedene Schemata für die Kodierung von Qualitätswerten verwendet, aber die Gemeinschaft hat sich auf das Format Sanger/Illumina 1.8+ (Phred+33) geeinigt, das hier beschrieben wird.

Phred-Score

Der Phred-Score (d.h., Punktzahl) ist definiert als

\[ Q = -10 \log_{10} \]

wobei \(p\) die Wahrscheinlichkeit ist, dass der entsprechende Basisaufruf falsch ist und \(Q\) die Phred-Punktzahl ist (gerundet auf den nächstliegenden ganzzahligen Wert). Der Phred-Score ist also eine einfache Transformation der Fehlerwahrscheinlichkeit, die eine einfache, aber einigermaßen platzsparende Kodierung darstellt (Tabelle 1).

Anmerkung: Phred-Basisqualitätswerte und ihre zugehörige Fehlerwahrscheinlichkeit (p) und Basisgenauigkeit. Eine Auswahl von Werten vom niedrigsten (0) bis zum höchsten (93) Phred-Score, der in einer FASTQ-Datei darstellbar ist, wird gezeigt.

ASCII-Kodierung von Phred-Scores

Die auf der letzten Seite beschriebene Transformation wandelt eine Wahrscheinlichkeit in einen ganzzahligen Wert zwischen 0 und 93 um. Die Werte werden in der FASTQ-Datei nicht als ein- oder zweistellige Zahl gespeichert, sondern als einzelnes Zeichen (char>), was wiederum zu einer erheblichen Verringerung der Dateigröße führt. Nichtsdestotrotz können Phred-Qualitäten von 0 bis 93 ein sehr breites Spektrum an Fehlerwahrscheinlichkeiten repräsentieren, das von 1,0 (100% Fehlerwahrscheinlichkeit oder einfach eine falsche Base) bis zu \(10^{-9.3}\) reicht, was einem extrem genauen Base-Call entspricht.

Um die Phred-Bewertungen als Zeichen zu speichern, werden die Scores in ASCII-Zeichen umgewandelt. ASCII (American Standard Code for Information Interchange) ist ein früher Zeichencodierungsstandard für die Darstellung von Zeichen in Computern und anderen Geräten, der erstmals 1963 veröffentlicht wurde. Die ASCII-Codes 0 bis 31 sind nicht druckbar. Der ASCII-Code 007 beispielsweise entspricht einem Steuercode, der ursprünglich gesendet wurde, um bei älteren Systemen eine elektromechanische Glocke zu läuten oder bei einigen Computern einen Systemwarnton abzuspielen. Das erste druckbare Zeichen ohne Leerzeichen ist der ASCII-Code 33, und der letzte druckbare ASCII-Code ist 126. FASTQ-Dateien verwenden also die ASCII-Codes 33-126 zur Codierung der Phred-Qualitäten von 0 bis 93 (Tabelle 2).

Illumina FASTQ-Dateibenennungsschema

Illumina verwendet ein Standardbenennungsschema für FASTQ-Dateien. Es ist nützlich zu verstehen, wie dieses Schema aufgebaut ist. Das allgemeine Schema für solche Dateien lautet

d.h., english

Schauen wir uns die folgende GIAB-Datei (Genome in a bottle) an:

NIST7035_TAAGGCGA_L001_R1_001.fastq.gz

Die Bestandteile dieses Namens sind:

1. sample_name: NIST7035. Dies ist der Probenname, der im Probenblatt für den Sequenzierungslauf angegeben ist.
2. barcode_sequence: TAAGGCGA. Dies ist die Nukleotidsequenz des molekularen Barcodes, der zur Markierung der Probe für das Multiplexing verwendet wird.
3. lane: 001. Die Lane-Nummer (1--8).
4. Spurnummer (d.h., read number): 1. Die Read-Nummer für Paired-End-Reads. R1 bedeutet Read 1, und für einen Paired-End-Sequenzierungslauf gibt es eine zusätzliche Datei mit R2 (Read 2), deren Name ansonsten genau dem Dateinamen für Read 1 entspricht
5. set_number: 001. Die maximale Dateigröße von FASTQ-Dateien wird mit der Befehlszeilenoption --fastq-cluster-count des Skripts configureBclToFastq.pl festgelegt, das Teil der Illumina \index{CASAVA} CASAVA-Software-Suite gehört. Wenn mehr Daten vorhanden sind, werden die Daten in separate FASTQ-Dateien mit der entsprechenden Dateigröße aufgeteilt (Um nur eine einzige FASTQ-Datei zu erstellen, kann eine ``0'' angegeben werden). Die verschiedenen Dateien, die derselben Probe/demselben Barcode/derselben Spur entsprechen, werden durch die mit 0 gefüllte dreistellige Set-Nummer unterschieden.

Bestimmte Illumina-Sequenzer verwenden andere FASTQ-Dateibenennungsschemata. Einzelheiten finden Sie in der Illumina-Dokumentation.

Paired-End-Sequenzierung

Beachten Sie, dass bei Paired-End-Läufen die übereinstimmenden FASTQ-Dateien genau die gleiche Anzahl von Reads aufweisen müssen und die Reads in beiden Dateien die gleiche Reihenfolge haben müssen. Dies kann mit den UNIX-Befehlen zcat und wc überprüft werden. Der Befehl cat liest Daten aus Textdateien und gibt deren Inhalt auf der Befehlszeilenschnittstelle aus, und der Befehl \verb+zcat+ tut dasselbe mit gzip-komprimierten Dateien. Der Befehl wc zählt Wörter, Zeilen und Zeichen in Textdateien. Wenn wir die Befehle wie folgt kombinieren, sehen wir, dass jede der beiden heruntergeladenen Dateien die gleiche Anzahl von Zeilen hat.

$ zcat NIST7035_TAAGGCGA_L001_R1_001.fastq.gz | wc -l
80812008
$ zcat NIST7035_TAAGGCGA_L001_R2_001.fastq.gz | wc -l
80812008

Wir können nun die Gesamtzahl der Zeilen durch vier teilen, um die Gesamtzahl der Reads zu erhalten (da jeder Read insgesamt vier Zeilen in der FASTQ-Datei einnimmt). Beachten Sie, dass wir nicht einfach nach Zeilen suchen können, die mit @ beginnen, da das ASCII-Symbol, das einem Phred-Score von 31 entspricht, ebenfalls @ ist und somit auch Qualitätszeilen mit diesem Zeichen beginnen können.

Read-Benenunng

Betrachten Sie die folgende Kennzeichnung, die von einem Read aus dem in Kapitel~\ref{ch:data} beschriebenen Exom-Datensatz NA12878 stammt.

@HWI-D00119:50:H7AP8ADXX:1:1101:2100:2202 1:N:0:TAAGGCGA

Die Syntax der Namenszeilen entspricht dem folgenden Schema:

@{instrument}:{rujn}:{flowcell_ID}:{lane}:{tile}:{x-pos}:{y-pos}\
    {read}:{is_filtered}:{control_number}:{index} 

Die im obigen Etikett gespeicherten Informationen sind in Tabelle 3 zusammengefasst. Der erste Teil dieser Bezeichnung, bis zum Leerzeichen, wird als Lesename oder Bezeichner verwendet.

Hinweis: Jedes Leseetikett speichert Informationen in einem Standardschema (Das Schema wurde mit CASAVA 1.8 und später verwendet). Wir haben die (optionale) UMI-Feld (Unique Molecular Identifier) weggelassen, da es nicht für die Exom- oder Genomsequenzierung verwendet wird.

S. auch Cook, 2010.

FASTQ-Analyse mit Rust

In diesem Abschnitt wollen wir eine einfache Rust-Version des Tools FastQC schreiben.

GitHub

GitHub ist eine Online-Plattform, die Entwicklern und Teams Werkzeuge zur Versionskontrolle und Zusammenarbeit an Softwareprojekten bietet. Sie basiert auf Git, einem weit verbreiteten Versionskontrollsystem, und erweitert dessen Funktionen um eine benutzerfreundliche Weboberfläche sowie zahlreiche zusätzliche Funktionen.

Falls Sie noch nicht mit GitHub gearbeitet haben, arbeiten Sie diesen Tutorial durch.

Praktikumsteilnehmer werden kleine Teams bilden, und am Entwicklen der Software zusammenarbeiten. Hierbei werden wir main'' und develop'' Branches verwenden. Teilnehmer sollen Feature-Branches kreieren und Pull-Requests erstellen, welche von anderen Teilnehmern zu kommentieren und ggf. anzunehmen sind.

Testgetriebene Entwicklung (Test-Driven Development; TDD)

Die TDD ist eine Methode zum Schreiben von Code, bei der ein automatisierter Testfall auf Unit-Ebene geschrieben wird, der fehlschlägt, dann wird gerade genug Code geschrieben, um den Test zu bestehen, dann werden sowohl der Testcode als auch der Produktionscode überarbeitet und dann mit einem anderen neuen Testfall wiederholt.

TDD in Rust

Wir werden das Erstellen von Unit-Tests in Rust während des Praktikums erläutern und demonstrieren. Es folgen Erklärungen zu einem einfachen Beispiel.

Übung 1

Die Definition des Phred-Scores in Bezug auf die Fehlerwahrscheinlichkeit wurde in diesem Abschnitt erklärt. Analog dazu lässt sich die Fehlerwahrscheinlichkeit aus der Phred-Punktzahl wie folgt berechnen.

\[ p = 10^{-Q/10}
\]

Konvertieren Sie die ASCII-Symbole für die Basen des folgenden Reads in Phred-Scores und berechnen Sie den durchschnittlichen Phred-Score dieser Basen. Verwenden Sie dann die Gleichung, um die durchschnittliche Fehlerwahrscheinlichkeit zu berechnen.

Wir werden im Laufe des Arbeitens mit FASTQ-Dateien überlegen (müssen), welche Datenstrukturen am besten geeignet sind, um die erforderlichen Analysen korrekt und effizient durchzuführen. Für diese Übung wollen wir damit anfangen, einige Funktionen zu schreiben, die den Qualitätsstring bewerten. Es folgt ein weiteres Bespiel.

@My-Illu:6:73:941:1973#0/1
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

Der folgende einfache Algorithmus kann verwendet werden, um den durchschnittlichen Qualitätswert zu berechnen.

Algorithm read_qual(qual_string)
    let n = len(qual_string);
    let qual_sum = 0;
    for qual in qual_string {
        let q = calculate_phred(qual)
        qual_sum += q
    } 
    return qual_sum / n

Wir richten die TDD für die function ``calculate_phred'' ein. Wir kreieren eine Datei namens fastq.rs und schreiben die folgende Funktion (die natürlich nur für ein bestimmtes Argument das richtige Ergebnis liefert!).

#![allow(unused)]
fn main() {
pub fn calculate_phred(qual: char) -> usigg {
    42
}
}

Wir richten dann eine Testfunktion ein (in der Regel lohnt es sich, das Testmodul am Ende der jeweiligen Code-Datei unterzubringen)

#![allow(unused)]
fn main() {
#[cfg(test)]

mod test {
    fn test_calculate_phred {
        let qual: char = '&';
        let expected:usign =	5; // Phred-Score für das Zeichen '&'
        let res = calculate_phred(qual);
        assert_eq!(expected, res);
    }

    // Alternativ können viele Ergebnisse auf einmal getestet werden:
    fn test_calculate_phred {
        let tests: Vec<char, usign> = vec![
            ('&', 5),
            ('+', 10)
        ];
        for test in tests {
            let res = calculate_phred(test.0);
            assert_eq!(test.1, res);
        }
    }
    // ggf. andere Testfunktionen
}
}

Wenn wir diesen Test zum ersten Mal durchführen, wird er natürlich einen Fehler melden. Wir arbeiten so lange an der Funktion ``calculate_phred'' bis die korrekt funktioniert. Wir arbeiten ggf. an der Funktion solange weiter, bis sie effizient ist. Dann gehen wir zur nächsten Funktion...Mit jeder neuen Funktion schreiben wir zunächst einen Test und dann die Funktion. Wir lassen jedes Mall alle Tests laufen. Das volle Konzept der TDD umfasst vieles mehr, aber für dieses Praktikum, wollen wir die hier skizzierte, einfache Vorgehensweise für alle Übungen anwenden.

Schreiben Sie eine korrekte Version der Funktion. Testen alle relevanten Argumente (kann man das "schlau" machen?).

Schreiben Sie dann eine Funktion für read_qual und testen Sie diese Funktion.

Übung 2

Schreiben Sie ein kleines Rust-Programm, um die in dem folgenden Datei-Namen enthaltenen Information zu erklären.

@HWI-D00107:50:H6BP8ACWV:5:2204:10131:51624 2:N:0:AGGCAGAA

Übung 3

Schreiben Sie ein kleines Rust-Programm, um die in dem folgenden Read-Namen enthaltenen Information zu erklären.

@Machine42:1:FC7:7:19:4229:1044 1:N:0:TTAGGC

\subsubsection*{Exercise 4} What does the quality character \verb+>+ mean? What is the associated error probability?

Übung 4

Stellen wir uns einen Read zu 1000bb vor, bei dem alle Basen einen Phred-Score von & aufweisen. Berechen Sie die erwartete Anzahl von Fehlern für den Read (unter der Unabhängigkeitsannehma).

FasQC in Rust

FastQC ist eine Applikation, die unter anderem eine Qualitätsbewertung von FASTQ-Dateien ausführt und die Ergebnisse in tabellarischer Form bzw. graphisch darstellt. Inzwischen gibt es viele neuere Applikationen, aber FastQC wurde früher von sehr vielen Gruppen sehr oft verwendet.

Zur Kontrolle Ihrer Ergebnisse können Sie FastQC hier herunterladen.

Ziel dieses Abschnitts ist es, die folgenden Funktionen zu als Rust-Applikation zu programmieren.

Befehlszeilenargumentenparser

Ihre Applikation soll den Pfad zu einer (Single-End) oder zwei (Paired-End) FASTQ-Dateien also Kommandozeilenargumente einlesen.

Wir empfehlen clap.

Beschreibende Statistiken

Als erstes überprüfen Sie ob die Datei existiert. Falls nicht, geben Sie eine verständliche Fehlermeldung aus.

Geben Sie die Informationen über den Dateinamen aus (s. hierzu letzte Seite)

Sequenzqualität (pro Base)

Geben Sie eine Tabelle aus mit der durchschnittlichen Qualität (Phred-Wert) für jede Position der Reads.

Sequenzqualität (pro Sequenz)

Geben Sie die durchschnittliche Qualität aller Reads aus.

Sequenzidentität (pro Base)

Geben Sie die durchschnittlichen Anteile von A, C, G, T für jede Position der Reads aus.

GC-Gehalt (pro Base)

Der GC-Gehalt definiert sich als prozentualer Anteil der Nukleotidbasen Guanin und Cytosin an der Gesamtheit der vier Basen (Adenin, Thymin, Guanin, Cytosin) der DNA. Sie den Wikipedia-Eintrag. Geben Sie den durchschnittlichen GC-Gehalt für jede Position der Reads aus.

GC-Gehalt (pro Sequenz)

Geben Sie den durchschnittlichen GC-Gehalt aller Reads aus.

Längenverteilung

Geben Sie eine Tabelle mit der Verteilung der Längen der einzelen Read aus (Anzahl und Prozentsatz).

SAM/BAM

Das SAM-Format (Sequence Alignment/Map) ist ein tabulatorgetrenntes Textformat, das ein universelles Format für die Speicherung von Alignments von NGS-Reads zu einem Referenzgenom sein soll.¹ BAM (Binary Alignment/Map) ist die komprimierte (binäre) Version von SAM. Das SAM/BAM-Format kann Alignments von allen wichtigen Alignern wie Bowtie 2 oder BWA speichern, unterstützt die Indexierung für eine schnelle Suche und Betrachtung mit dem IGV und wird von Software zur Variantenerkennung weitgehend unterstützt.

Es gibt eine umfangreiche Online-Dokumentation des Formats, und wir werden hier nicht versuchen, alle Details zu behandeln. Stattdessen werden wir die verschiedenen Komponenten des Formats anhand von kleinen Beispielen intuitiv erklären. Darüber hinaus zeigen wir, wie man die Eigenschaften der durch SAM/BAM-Dateien dargestellten Alignments analysiert. Zu diesem Zweck werden wir hauptsächlich SAMtools verwenden.

SAM-Grundlagen

SAM-Dateien bestehen aus einem optionalen Header-Abschnitt, gefolgt von einem Alignment-Abschnitt. Die Kopfzeilen beginnen mit `@'. Im Alignment-Abschnitt der SAM-Datei stellt jede Zeile ein Alignment für einen Read dar und umfasst mindestens 11 Pflichtfelder (Tabelle 1).

Jede Zeile im Alignmentsabschnitt einer SAM-Datei umfasst 11 Pflichtfelder.

Das SAM/BAM-Format soll so allgemein wie möglich sein und verwendet das Wort template für DNA-Fragmente (englisch: fragments, inserts). Das SAM-Format verwendet das Wort segment, um sich auf eine zusammenhängende Sequenz oder Teilsequenz zu beziehen; zum Beispiel können die beiden Reads eines Lesepaares als zwei Segmente bezeichnet werden.

Das Feld QNAME speichert den Namen der Abfrage-(query)-Sequenz, in der Regel ein Read.¹ POS gibt die Position auf dem Chromosom an, an der der Read ausgerichtet wurde. Das Feld FLAG wird weiter unten erklärt.

Das Feld MAPQ stellt die Mappingqualität, die vom Aligner zugewiesen wurde, und spiegelt das Vertrauen wider, mit dem der Read auf die angegebene Position abgebildet werden konnte. Die CIGAR-Zeichenkette bietet eine kompakte kompakte Darstellung des Alignments, die im Folgenden erläutert wird.

Die Felder RNEXT, PNEXT und TLEN werden nur für Paired-End-Reads verwendet. Daher wird bei Single-End Reads RNEXT auf ``*'' (für not-applicable) gesetzt. PNEXT und TLEN werden auf 0 gesetzt. Schließlich enthält das Feld SEQ die Nukleotidsequenz des Reads und das Feld QUAL zeigt die Basenqualitätswerte für jede Position an.

Paired-End-Sequenzen

In Abbildung 2 werden zwei Paare bestehend aus vier Reads der Länge 10 werden auf eine Referenzsequenz der Länge 40 abgebildet. Beachten Sie, dass Mitglieder desselben Paars haben denselben QNAME und werden auf verschiedene Stränge abgebildet. Read 1 wird dem dem Vorwärtsstrang an Position 2 zugeordnet und bildet ein Paar mit dem Read, der dem Rückwärtsstrang an Position Strang an Position 27 zugeordnet ist. RNEXT wird auf = gesetzt, da der zweite Read des Paares des Paares auf dieselbe Referenz, d. h. chrE, abgebildet wird, und PNEXT wird auf 27 gesetzt, weil der Reverse-Read auf diese Position abgebildet wird (der Reverse-Read ist an den Positionen 27-36 ausgerichtet). Daher wird TLEN auf 35 gesetzt, denn dies ist der Abstand zwischen der äußersten linken und der äußersten rechten der am weitesten rechts gelegenen Base der beiden Reads. Der Eintrag TLEN für den Reverse-Read wird mit einem Minuszeichen dargestellt.

Für Paired-End-Reads werden andere bitweise Flags verwendet als für Single-End-Reads. Die Flags 99 und 147 zeigen an, dass beide Mitglieder des Paares korrekt auf den Vorwärts- und den Rückwärtsstrang gemappt wurden.

CIGAR

Eine CIGAR (Concise Idiosyncratic Gapped Alignment Report) String besteht aus einer Reihe von Operationslängen sowie den Operationen, die beschreiben, wie genau ein Read an die Referenzsequenz angeglichen wurde (Abbildung 1).

Die drei wichtigsten CIGAR-Operationen sind M (Match/Mismatch), I (Insertion) und D (Deletion), aber aber es gibt noch sechs weitere Operationen, von denen einige im Folgenden erläutert werden. Alle CIGAR Operationen sind in Tabelle 1 zusammengefasst. Eine CIGAR-Zeichenkette kann man sich als eine Reihe von Operationslängen plus die CIGAR-Operation, die eine abstrakte und kompakte Darstellung beliebiger Ausrichtungen ermöglicht. Betrachten Sie beispielsweise das folgende Alignment, in dem die Abfragesequenz (Q) an die Referenzsequenz (R) angeglichen wird, beginnend an Position 7 der Referenzsequenz.

AGCATGTTAGATAA--GATAGCTGG R
      || ||||| |||| ||||
------TTGGATAAAGGATA-CTGG Q

Es gibt 8 anfängliche Übereinstimmungen, dann eine Einfügung von 2 Nukleotiden in die Abfrage, dann vier weitere Übereinstimmungen, dann eine Löschung von einer Base in der Abfrage und dann vier weitere Übereinstimmungen. Wir können dies als 8M2I4M1D4M darstellen.

Man beachte, dass die dritte Base der Abfragesequenz (G) sich tatsächlich von der Base der Referenzsequenz (A) unterscheidet, an die sie angeglichen wurde. Nichtsdestotrotz handelt es sich um eine "Alignment match" (Ausrichtungsübereinstimmung) (M), obwohl es keine Sequenzübereinstimmung ist.

Übung 1

Nehmen Sie eine Referenz- und eine Querysequenz sowie den CIGAR-String, der das Alignment spezifiziert (z.B. von Abb. 1). In unserem Beispiel wären die drei Argumente also

• Query (read_4): GCATCTGTTG
• Ref: GCATACTGTTG
• CIGAR: 4M1D6M

Das Ziel ist es, eine Funktion zu schreiben, welche das entsprechende Alignment ausgibt, also in unserem Beispiel

GCAT-CTGTTG
|||| ||||||
GCATACTGTTG

Übung 2

Probieren Sie Ihre Funktion mit "realistischen" BAM-Dateien -- wir haben eine TODO -- "kleine" BAM-Datei erstellen.

Clipped Reads

Wichtige Mapping-Tools wie BWA versuchen, einen Teil eines Reads zu kartieren, wenn sie den Read nicht in voller Länge auf das Referenzgenom abbilden können. In diesen Fällen kann der nicht zuzuordnende Teil durch einen Prozess namens Clipping aus dem Alignment ausgeschlossen werden.

Mit soft-clipping, das in der CIGAR-Zeichenkette mit einem >S angezeigt wird, werden die abgeschnittenen Sequenzbasen nicht aus der SEQ-Zeichenkette entfernt, aber nicht von Varianten-Callern verwendet und nicht in Viewern wie IGV (Integrativer Genomik-Viewer)angezeigt.

Das Hard-Clipping (H) ähnelt dem Soft-Clipping (S), unterscheidet sich aber dadurch, dass die hart geclippte Teilsequenz nicht im Alignment-Datensatz enthalten ist (Abbildung 1).¹

Ref:   GTTCCTAGGAACAGCACAATTTCTAGATACAATCAT
Read1:    CCTAGGAACAGCACAATTTCTAGATACAATCAT
Read2:     ggtcacatgattgtATTTCTAGATACAATCAT

Das entsprechende Alignment von Read1 und Read2 mit Ref. Die soft-clipped (nicht ausrichtbaren) Basen von Read2 sind in Kleinbuchstaben dargestellt.

Mapping-Qualität

Das fünfte Feld einer SAM-Datei gibt die Mapping-Qualität des Read an. Der MAPQ-Wert (MAPping Quality) spiegelt die Wahrscheinlichkeit wider, dass der Read an der falschen Position im Genom ausgerichtet ist. Wie der Wert für die Basenqualität handelt es sich dabei um eine Phred-skalierte posteriore Wahrscheinlichkeit, dass die vom Aligner angegebene Mapping-Position falsch ist

\[ MAPQ = -10\log_{10} P(\mathrm{mapping position wrong) \]

Die MAPQ wird auf die nächste ganze Zahl gerundet. Beispielsweise würde einer von tausend Reads mit einem MAPQ von 30 als falsch ausgerichtet vorhergesagt. Ein MAPQ-Wert von 255 bedeutet, dass die Zuordnungsqualität nicht verfügbar ist.

Eine der größten Herausforderungen für Read-Alignment-Algorithmen ist das Vorhandensein von sich wiederholenden Sequenzen im menschlichen Genom (repeats), so dass ein kurzer Read an zwei oder mehr Positionen gleich gut ausgerichtet sein kann. In diesem Fall weist BWA dem Read eine MAPQ von Null zu und wählt eine der Positionen nach dem Zufallsprinzip aus.

Man sollte sich der Tatsache bewusst sein, dass verschiedene Aligner unterschiedliche Methoden zur Berechnung der Mapping-Qualität verwenden und die von verschiedenen Alignern erzielten MAPQ-Scores im Allgemeinen nicht direkt vergleichbar sind.

Übung 1

Berechnen Sie die durchschnittliche Mapping-Qualität aller Reads in einer BAM-Datei. Berechnen Sie den Anteil der Reads mit einer Mindestqualität von Q30.

Das SAM-Bitflag

Das zweite Feld jedes Datensatzes in einer SAM-Datei stellt ein Bitfeld mit Werten aus 12 Bitflags~(Tabelle 1) dar.

Das SAM-Format verwendet 12 Bitflags, von denen jedes den Wert 1 (Ja, Wahr) oder 0 (Nein, Falsch) haben kann. Die Bitflags können unabhängig voneinander kombiniert werden und werden in SAM-Dateien als entsprechender dezimaler Ganzzahlwert angezeigt.

Ein Bitfeld wird verwendet, um eine Reihe boolescher Werte (Ja/Nein) kompakt zu speichern. In unserem Fall möchten wir 12 Ja/Nein-Attribute zu jedem Lesevorgang speichern. Theoretisch könnten wir zwölf separate char-Werte (jeweils ein Byte) speichern, aber wenn wir beachten, dass die 12 Werte als einzelne Bits eines Bitfelds gespeichert werden können, können wir eine erhebliche Menge an Speicherplatz einsparen. Tabelle 1 zeigt die zwölf Bitflags und Abb. 1 zeigt ein Beispiel dafür, wie die Bitflags als dezimale Ganzzahlwerte dargestellt werden.

Oft ist es zweckmäßig, die hexadezimale Schreibweise („hex“) zu verwenden, bei der die Basis nicht 10 oder 2, sondern 16 ist; Die Ziffern werden durch 0–9 dargestellt, gefolgt von A, B, C, D, E und F. Um mit SAM-Dateien völlig vertraut zu sein und wenn Sie Programme wie SAMtools verwenden, um damit zu arbeiten, ist es hilfreich zu verstehen, wie man konvertiert zwischen diesen drei Zahlendarstellungen.

Ein Wert von 99 zeigt zum Beispiel folgende Ergebnisse an:

• der Read mehrere Segmente in der Sequenzierung hat,¹ d.h., 0x1= \(1\times 16^0=1\)
• jedes Segment des Reads konnte richtig gemappt werden (korrektes Mapping beider Reads des Read-Paares), d.h., 0x2=\(2 \times 16^0=2\)
• der Mate-Read auf den Rückwärtsstrang gemappt ist, d.h., 0x20=\(2\times16^1=32|))
• der Read, auf den im aktuellen SAM-Datensatz Bezug genommen wird, ist das erste Segment im Tempalte (der erste Read im dem Paar, d.h., 0x40=\(4\times 16^1=64\))

Um zu verstehen, wie die Bitfelder dargestellt werden, ist es hilfreich, sich an die Position zu erinnern Die Notation von Zahlen funktioniert. In der bekannten Dezimalschreibweise steht die Zahl 453 für \(4\times 10^2 + 5\times 10^1+3\times10^0= 400 + 50 + 3 = 453\). In binärer Schreibweise steht eine Zahl wie 101101 für \(1\times 2^5+0\times 2^4+1\times2^3+1\times 2^2+0\times 1^1+1\times 2^0 = 32 + 8+4+1= 45\) - vgl. hierzu Abbildung 1 und Tabelle 1.

Bit-Operationen können verwendet werden, um Lesevorgänge in einer SAM-Datei zu filtern. Das SAM-Format verlangt zum Beispiel, dass für jeden Lesevorgang eine und nur eine der zugehörigen Zeilen folgende Bedingungen erfüllt

FLAG & 0x900 == 0

Diese Zeile wird als primäre Zeile (primary line) des Reads bezeichnet. Der hexadezimale Wert 0x900 ist gleich 0x800 plus 0x100, was bedeutet, dass die Bit-Flags für das supplementäre Alignment (0x800) und das sekundäre Alignment (0x100) gesetzt sind. Erinnern Sie sich, dass & für die bitweise logische UND-Verknüpfung eines Bitpaares steht. Wenn beide Bits 1 sind, ist das Ergebnis der Operation 1, ansonsten 0. Zum Beispiel

      1010       
AND   1100       
     =1000      

Wenn also ein Read entweder als sekundäres Alignment² oder ein \index{supplementary alignment} supplementary alignment,\footnote{Nach der SAM-Spezifikation ist ein chimäres Alignment ein Alignment eines Reads, das nicht als lineares Alignment dargestellt werden kann. In der Regel bedeutet dies, dass ein Read aus mehreren Segmenten besteht, die sich an verschiedenen Teilen des Genoms ausrichten. Die Segmente haben keine großen Überlappungen. Eines der Segmente wird als das repräsentative Alignment betrachtet, die anderen werden als ergänzend bezeichnet und haben das Flag für ergänzendes Alignment gesetzt. dann erhalten wir FLAG & 0x900 != 0.³ Das SAM-Format besagt, dass eine beliebige Anzahl von Zeilen sekundäre oder ergänzende Ausrichtungen darstellen kann, aber dass nur eine Zeile die primäre (repräsentative) Ausrichtung sein kann.

samtools

SAMtools kann Alignments filtern, um nur die alignierten Reads auszugeben, die bestimmten Kriterien entsprechen, einschließlich der Überlappung mit Regionen, die in einer BED-Datei (-L) angegeben sind, Zugehörigkeit zu bestimmten Lesegruppen (-r, -R), Herkunft aus einer bestimmten Bibliothek oder eine bestimmte Mindestanzahl von CIGAR-Basen haben, die die Abfragesequenz verbrauchen.¹

Die Option -f gibt nur Alignments aus, bei denen alle Bits in INT im Feld FLAG gesetzt sind (INT kann als Dezimalzahl oder als hexadezimale Zahl angegeben werden, die mit '0x' beginnt). Dies führt eine UND-Verknüpfung mit allen Bits durch, die in INT gesetzt sind, was bedeutet, dass Sequenzen nur angezeigt werden, wenn mindestens die entsprechenden Bit Positionen in FLAG auf 1 gesetzt sind. Die Option -F gibt keine Alignments aus, bei denen Bits, die in INT gesetzt sind, im Feld FLAG. Die Option -F überspringt daher alle Reads, für die FLAG & INT != 0 ist. Um also nur primäre Reads auszugeben, können wir den folgenden Befehl verwenden.

$ samtools view -F 0x900 NIST7035_aln.bam

Die mit diesem Befehl ausgegebenen Reads enthielten insgesamt 38 verschiedene Werte für das Bitfeld, von denen keiner die Bits 0x100 oder 0x800 enthielt. Mit dem folgenden Befehl werden alle sekundären Alignment-Reads extrahiert.

$ samtools view -f 0x100 NIST7035_aln.bam

Die Funktion flagstat von SAMtools bietet eine Zusammenfassung der Anzahl der Datensätze, die jedem der Bit-Flags von Tabelle 1 entsprechen. Die Auswertung unseres Beispieldatensatzes zeigt beispielsweise Folgendes.

$ samtools flagstat NIST7035_aln37.bam 
35210329 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
35210007 + 0 mapped (100.00%:-nan%)
35210329 + 0 paired in sequencing
17605154 + 0 read1
17605175 + 0 read2
34991352 + 0 properly paired (99.38%:-nan%)
35209746 + 0 with itself and mate mapped
261 + 0 singletons (0.00%:-nan%)
19085 + 0 with mate mapped to a different chr
9718 + 0 with mate mapped to a different chr (mapQ>=5)

SAM-TAGs: Optionale Felder

Auf die 11 Pflichtfelder können beliebig viele optionale Felder folgen. Alle optionalen Felder folgen dem TAG:TYPE:VALUE-Format, wobei TAG eine zweistellige Zeichenfolge, TYPE einer von sechs Datentypen (Tabelle 1) und VALUE der tatsächliche Wert ist.

Verschiedene Aligner verwenden unterschiedliche optionale Felder. Leser finden eine vollständige Dokumentation allgemeiner Felder bei SAMtools Website¹ und Aligner wie BWA stellen ebenfalls Dokumentation zu den von ihnen verwendeten Feldern bereit. Hier stellen wir einige ausgewählte Bereiche vor, um den Lesern einen Eindruck davon zu vermitteln, wie diese Bereiche in der Praxis eingesetzt werden.

Das NM-Feld