# Vorlesung, 8. Woche
###### tags: `Vl Datenanalyse 21/22`
## Hochdurchsatz-Sequenzierung
(High-throughput sequencing, HTS; Next-generation sequencing, NGS)
Illumina-Sequenzier-Gerät:
Moderne Geräte können mehrere 100 Millionen DNA-Fragmente pro Tag sequenzieren.
Vorgehen:
Library Prep:
- Eine Probe mit sehr vielen verschiedenen DNA-Fragmenten (mit max. ~500 bp Länge) wird vorbereitet.
- Sequenzier-Adapter werden an die Fragmente anligiert (oder durch Template-Switch-PCR o.ä. angefügt).
- Durchführung einer PCR, um von jedem Fragment viele Kopien zu erhalten. Man erhält die sog. "Library".
Cluster generation:
- Die Library wird in die Flowcell gefügt; die DNA haften sich am Primer-Rasen auf der Oberfläche an.
- Die Bridge-PCR fürt dazu, dass jedes DNA-Molekül zu einem kleinen "Cluster" aus Kopien wächst.
Sequencing-by-Synthesis:
- Im Sequenzierer wird die Sequenz der Cluster abgelesen, indem Schritt für Schritt floureszent markierte Nukleinsäuren angehängt werden.
- Aus automatisch erstellen Mikrographien wird die Sequenz abgelesen.
## RNA-Sequencing: Überblick
Soll RNA statt DNA sequenzietrt werden, wird zunächst eine Reverse Transkription (RT) durchgeführt.
Wikipedia:
## FASTQ-Dateien
Der Sequenzierer spuckt Textdateien aus, die so aussehen:
```
@ERR649019.76 76/1
TTGCACCTCGCCCATTTGGAGATATTAGTGGAGAATTGGTCTGGTCTTGATCTCTGGGTAGAGATGTGAGTGGGAGTTATCAATTTGGGAGTCATCCGAA
+
CCCFFFFFHHHHHJJJJJJGIHIJJJJJGHIHGIJJJJIFGIJJIIIJJIJJJJIIJJHHJJJJJJGIJJEHHHHD?CFEEEEDEEEDDDDBDDED>@DB
@ERR649019.77 77/1
TGTGCCTGTAGTCCCAGCTACTCGGGAGGCTGAGGTGGGAGGATCGCTTGAGCCCAGGAGTTCTGGGCTGTAGTGCGCTATGCCAGATCGGAAGAGCACA
+
@C@FFFFFFHHFFHIIHIJEIHIJJJHIIIJIIIJ9BFFFHIJIJJJIIII@CHHEEDFF?CEEEEEDDDDDDDDCDDDDDDDDDDDDDDDDDDDDCDDC
@ERR649019.78 78/1
NGGTAGATCGGAAGAGCACACGTCTGAACTCCAGTCACGGCTACNNNNNNNATGCCGNCTTCTGCTTGAAAAAANACACAAATATANNNNNNNTNGNNAG
+
#1:BDDDDDDDDDIAEEEEF8@FEBBEEEIEIEBB?BD?A@6?B########################################################
```
Jeder Read wird durch vier Zeilen beschrieben:
- `@`, gefolgt von der Read-ID (selten relevant, enthält oft die Seriennummer des Sequenzierers und der Flowcell, und die Position auf der Flowcell)
- der eigentliche Read
- ein einsames `+`
- der Quality-String, der kodiert, wie sicher der Sequenzierer war, dass die "Base calls" (also die Buchstaben der Sequenz) richtig sind.
Eine typische FASTQ-Datei enthält viele Millionen von Zeilen.
Beim Paired-End-Sequenzieren macht der Sequenzierer zwei Durchläufe und prodziert zwei FASTQ-Dateien, die die Cluster in derselben Reihenfolge beschreiben. Der jeweils erste Read ("R1") beschreibt das 3'-Ende des mRNA-Fragments, der zweite Read ("R2") das 5'-Ende (oder, bei manchen Library-Prep-Verfahren, auch umgekehrt).
### FASTQ-Dateien der OSCC-Studie
Im [Paper zur OSCC-Studie](https://doi.org/10.18632/oncotarget.5529) steht am Ende: "Our sequencing data has been made publicly accessible via the European Nucleotide Archive under study accession PRJEB7455 (http://www.ebi.ac.uk/ena/data/view/PRJEB7455)."
Dort finden wir für jede der 57 Proben je zwei FASTQ-Dateien (Paired-End), und eine Tabelle, die die Run-IDs (Dateinamen der FASTQ-Dateien) den Proben-IDs zuordnet.
## Read Alignment
Für jedes Read-Paar möchten wir nun die Stelle im Genom finden, von dem das Fragment stammt.
Für einzelne Reads kann man BLAST verwenden. Da wir aber Millionen von Reads haben, verwenden wir einen speziellen "Short-Read-Aligner".
Die meisten Aligner startet man von der Kommandozeile; wir verwenden hier einen, den man von R aus starten kann: den "Subread"-Aligner, entwickelt von Liao, Smyth und Shi.
Rsubread ist ein Teil von Bioconductor, einer Sammlung von R-Paketen für Bioinformatik.
Daher installieren wir zunächst den Bioconductor-Paket-Manager:
```
install.packages( "BiocManager" )
```
Mit diesem können wir nun andere Bioconductor-Pakete installieren (die man nämlich nicht direkt mit `install.packages` installieren kann):
```
Biocmanage::install( "Rsubread" )
```
### Genom-Sequenz
Wir laden uns von [Ensembl](www.ensembl.org) eine FASTA-Datei mit der aktuelle Fassung des menschlichen Genoms herunter: "GRCh38", d.h. Assembly-Version 38, veröffentlicht vom "Genome Reference Council" (GRC), für "human" (h). Wir verwenden die Datei `Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz`.
Die ist eine mit GZip (`.gz`) komprimierte FASTA-Datei (`.fa`).
Eine FASTA-Datei ist eine Text-Datei mit DNA-Sequenzen. Jede Sequenz, z.B. ein Chromosom, beginnt mit einer Zeile mit `>` und dem Sequenznamen (also z.B. `>1` oder `>chr1` für Chromosom 1), gefolgt von (oft vielen Tausend) Zeilen mit der DNA-Sequenz.
Die FASTA-Datei mit dem "primary assembly" enthält die Sequenzen der 22 Autosomen, der beiden Sex-Chromosomen, der Mitochondrien, sowie ...
### Index-Bau
Damit der Aligner die Reads schnell finden kann, muss er das Genom in einer übersichtliche Form bringen. Diese Vorbereirung nennt man "building a genome index"
Wir fordern RSubread auf, aus der FASTQ-Datei den Index zu bauen und unter dem Namen `GRCh38` zu speichern.
```
library( Rsubread )
buildindex( "GRCh38", "~/Downloads/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz",
gappedIndex=TRUE, indexSplit=FALSE, memory=4000 )
```
Das dauert etwa 30 Minuten und erzeugt einige Dateien, deren Namen mit `GRCh38` beginnt.
### Alignment
Nun können wir für eine Library ein Alignment erzeugen:
```
align( "GRCh38",
"ERR649018_1.fastq.gz", "ERR649018_2.fastq.gz",
output_file="ERR649018.bam" )
```
Wir haben angegeben: Den Namen des Index, den wir gebaut haben, den Namen der beiden Fastq-Dateien, die zum Paired-End-Run ERR649018 gehören, und den Namen, den die Datei mit den Alignments bekommen soll. Laut Tabelle ist Run ERR649018 die zu Probe "PG038-D" gehörige Library.
Das Alignment dauert knapp eine halbe Stunde.
Wir erhalten am Ende diese Zusammenfassung:
```
//================================ Summary =================================\\
|| ||
|| Total fragments : 32.177.535 ||
|| Mapped : 24.716.114 (76,8%) ||
|| Uniquely mapped : 23.060.524 ||
|| Multi-mapping : 1.655.590 ||
|| ||
|| Unmapped : 7.461.421 ||
|| ||
|| Properly paired : 14.014.830 ||
|| Not properly paired : 10.701.284 ||
|| Singleton : 1.701.060 ||
|| Chimeric : 36.301 ||
|| Unexpected strandness : 5.035 ||
|| Unexpected fragment length : 121.008 ||
|| Unexpected read order : 8.837.880 ||
|| ||
|| Indels : 271.812 ||
|| ||
|| Running time : 29,5 minutes ||
|| ||
\\============================================================================//
```
Subread bietet zwei Aligner an: `align` ist schneller, `subjunc` untersucht hingegen genau, wie ein Read aliniert werden muss, wenn er eine Splice-Junction enthält, also im Referenz-Genom über ein Intron springt.
### SAM-Dateien
Das Ergebnis von `align` ist eine sog. "binary sequence alignment/map"-Datei (BAM-Datei, oder binäre SAM-Datei). Mit einer speziellen Software (samtools) können wir die Datei ansehen, indem wir in der Kommandozeile (nicht in R) eingeben:
```
samtools view ERR649019.bam | less -S
```
SAM-Dateien sind eine Tabelle mit folgenden Spalten
- Read-ID des Reads, der aliniert wurde
- Flag-Feld (kodiert verschiedene Information, z.B. "Strand")
- Referenz-Sequenz-Name aus der FASTA-Datei, also z.B. Chromosom, auf dem die Read-Sequenz gefunden wurde
- Position, wo die Sequenz im Chromosom gefunden wurde
- Alignment-Qualität: sagt, wie sicher der Aligner ist
- CIGAR-String: kodiert, ob es Lücken o.ä. im Alignment gab, z.B. von Introns
- Informationen zum Alignment des Partner-Reads (anderes Ende des Fragments)
- Sequenz des Reads
- Quality-String des Reads
- Extra-Felder mit Zusatz-Informationen
### BAM-Dateien im Genom-Browser
Spezielle Programme wie z.B. der am Broad Institute entwickelte Integrated Genome Viewer (IGV) erlauben, SAM-Dateien grafisch anzusehen.
**TO DO**: Image
Damit IGV uns zeigen kann, wo die Gene auf den Chromosomen liegen, braucht es Genom-Annotationen.
## Genom-Annotation in GTF-Dateien
Das Genom ist die DNA-Sequenz. Die Genom-Annotation sagt und, wo die Gene liegen, d.h., wo Transkription zu RNA startet und endet, und wo Introns herausgespleißt werden.
Wir laden uns die von Ensembl für GRCh erstellte Annotation von der Ensembl-FTP-Seite herunter: Homo_sapiens.GRCh38.104.gtf.gz
Die Zahl "104" hier ist die Ensembl-Version: Ensembl revidiert die Annotationen etwa einmal pro Jahr und zählt dann diese Nummer hoch.
Eine GTF-Datei ist eine TSV-Datei. Wir können Sie also in R oder Excel ansehen.
Für jedes "Feature" gibt es eine Zeile. Features sind hierarchisch: Ein "gene" enthält mehrere Features vom Typ "transcript", diese sind aus Features vom Typ "exon" aufgebaut, und Exons enthalten Teile wie z.B. "CDS" (coding sequence), "three_prime_utr", und "five_prime_utr".
Für jedes Feature is Chromosom, Start- und End-Position (in Basenpaaren), und Strang angegeben, sowie ein Attribut-Feld, das Namen und ID des Gens und ggf. Transkripts enthält, Transkript-Biotype (coding, noncoding, rRNA, tRNA, etc. ) und einiges mehr.
Die Ensembl-Gen-ID beginnt beim Menschen mit "ENSG" und bezeichnet ein Gen auch dann eindeutig, wenn das Human Genome Naming Committee (HGNC) keinen verbindlichen Namen vergeben hat.
## Reads zählen
Nun möchten wir für jede Library ermitteln, wie viele Reads auf jedes Gen gefallen sind.
Dazu verwenden wir die Funtion `featureCopunt` aus dem `Rsubread`-Paket
```
featureCounts(
c( "ERR649018.bam", "ERR649019.bam" ),
annot.ext="Homo_sapiens.GRCh38.104.gtf",
isGTFAnnotationFile=TRUE, isPairedEnd=TRUE ) -> fc
```
Wir geben an:
- einen Vektor mit den Dateinamen aller BAM-Dateien (hier nur zwei, aber eigentlich alle 57)
- Die Annotation, hier als GTF-Datei
- die Information, dass es Paired-End-Libraries waren
Wir erhalten in `fc`
- die Count-Matrix `fc$count`, mit eienr Spalte für jede BAM-Datei und einer Zeile für jedes Gen aus der GTF-Datei
- in `fc$stat` Informationen darüber, wie viele Reads einem Gen zugeordnet werden konnten, und wie viele nicht
- in `fc$targets` und `fc$annotation` Informationen über unsere Eingabedaten
## Proben-Tabelle
Damit wir mit der Count-Matrix weiter arbeiten können, müssen wir eine Proben-Tabelle erstellen, die für jede Probe eine Zeile enthält mit den "Metadaten" der Probe: der Dateiname, die Patienten-Id, der Gewebe-Typ
## Differential gene expression analysis
In der Hausaufgabe haben wir nun "per Hand" mittels Tidyverse nach differentiell exprimierten Genen gesucht. Wir können das auch automatisiert machen, z.B. mit "limma-voom".
[...]
## Überblick behalten
[...]
Sign in with Wallet
Connect another wallet