Introductie

Dit project gaat vooral om het vinden van verschillen tussen een behandelconditie en een controleconditie op basis van whole-genome sequencing data. Daarvoor is het belangrijk om goed te kijken naar de kwaliteit van de ruwe reads, daarna de reads te aligneren tegen een referentiegenoom en vervolgens varianten te vinden en te interpreteren volgens de originele paper.

De gekozen dataset komt uit een Drosophila-experiment waarin meerdere condities met elkaar vergeleken worden. In dit geval is de Constant RNAi sample de tumor conditie en de Control sample de normale controle. De kern van de vraag is dus niet alleen: “zijn er varianten?”, maar ook: “welke afwijkingen zijn relevant en consistent met de experimentele opzet?” waarbij de varianten ook geannoteerd zullen moeten worden.

Centraal binnen het project staan de volgende voorwaarden: de data moet goed gekarakteriseerd zijn en de analysepijplijn moet reproduceerbaar en transparant worden uitgevoerd.

Ik werk samen met Vani Rembet, Dennis Kuiper en Kevin Matahelamual. Dit logboek dient als een document waarin de keuzes, problemen en oplossingen binnen het project terug te volgen zijn.

Project Layout

Alle logboeken, supplementaire code en visualisaties zijn te vinden op: Github (m.u.v. dit logboek)

De gebruikte data is te vinden op assemblix: /students/2026-2027/Thema05/GenTeam1/

Overige analyseresultaten die gevonden zijn zonder inbreng van de rest van het team te vinden onder hetzelfde pad + /andere_subsets/j_ss/.

Pipeline overview

Stap Tool / methode Doel Verwijzing binnen logboek
1 Artikel selectie Gepubliceerde dataset selecteren 07/09 - 15/09
2 Data download en QC SRA-download + steekproefcontrole 14/09 - 16/09
3 fastp kwaliteitscontrole en trimming 21/09 - 26/09
4 BWA-MEM aligneren tegen dm6 21/09 - 26/09
5 Picard / read groups / MarkDuplicates BAM-validatie en duplicaten markeren 21/09 - 23/09
6 Manta + Strelka2 somatische SNV/InDel-calling 26/09 - 05/10
7 CNVnator / CNVkit / Delly copy-number / structural variation 28/09 - 07/10
8 Resultaten en interpretatie variant- en CNV-bespreking 07/10 - 09/10

Gebruikte samples

Sample Run accession Condition Type Opmerking
Tumor / treated sample SRR22984524 Constant RNAi whole-genome WGS gebruikt als “tumor” in de somatische vergelijking
Control / normal sample SRR22984526 Control whole-genome WGS gebruikt als normale controle
- SRR22984519 Constant RNAi chip_antibody (active motif 39155) Foutief\(^{1}\)

\(^{1}\)Eerste experimenten met SRR22984519 zijn foutief, omdat die sample niet bij de juiste condition/whole-genome groep hoorde.

Onderzoeksvraag

De centrale onderzoeksvraag voor dit project is:

“Kunnen we met een reproduceerbare whole-genome analysepijplijn varianten of genomische afwijkingen detecteren die de behandelde toestand onderscheiden van de controleconditie?”

Met andere woorden: we willen bepalen of de data in de tumor conditie relevante genomische afwijkingen laat zien om een vergelijking met de controlestaat te kunnen maken. De keuze van de juiste samples en het volgen van de originele pipeline is hierbij essentieel.

Days

07/09

Vandaag begon het project. Ik wil eerst een geschikte dataset vinden, de repo opzetten en de basisstructuur van het project duidelijk krijgen zodat de rest van het project goed kan verlopen.

Plan:

  • Een geschikte dataset en artikel vinden voor genomics en transcriptomics.
  • De Git repo in orde maken.
  • Trello opzetten en de communicatiekanalen voorbereiden.

Gedaan:

  • Ik heb samen met de groep meerdere artikelen en datasets bekeken.

  • De volgende artikelen werden besproken en beoordeeld:

    • Afgekeurd artikel 1
      • Commentaar: hier leek geen standaard sequencing te zijn gebruikt, maar vooral Hi-C data. Dat is niet geschikt voor een genomics-variantanalyse. Het RNA-seq gedeelte was wel bruikbaar, maar niet als basis voor het genomics-gedeelte.
    • Afgekeurd artikel 2
      • Commentaar: de genomics-data was goed, maar de opzet was niet direct reproduceerbaar. Er was gebruik gemaakt van een uitgebreide GATK-workflow en Copy Number Variation, wat voor ons te complex was. RNA-seq was wel mogelijk, maar minder standaard.
    • Afgekeurd artikel 3
      • Commentaar: de data leek goed, maar het artikel ontbrak en de GEO-pagina gaf “Citation Missing”.
    • Goedgekeurd artikel 4 (Parreno et al. 2024)
      • Commentaar: dit leek het meest geschikt. Beide datasources waren aanwezig, de workflow was deels reproduceerbaar en zowel de genomics- als transcriptomicsdata waren beide goed te gebruiken. GATK was wel gebruikt, maar wij vonden dat we dit konden vervangen met alternatieve tools.
    • Mogelijk artikel 5
      • Commentaar: nog niet definitief beoordeeld, maar relevant als mogelijke alternatieve optie. \(\\\)

  • SSH-keys voor GitHub-toegang via het Bioinf-netwerk voor mij, Vani en Dennis ingesteld.

  • Eerste commit aan de gedeelde repo, inclusief dit logboek en een subfolder voor de logboeken v.d. anderen.

  • Het abstract van het goedgekeurde artikel is deels gelezen als eerste stap in het bepalen van de pipeline.

Conclusie: Op deze dag bekijk ik meerdere artikelen en datasets. Artikel 4 lijkt het meest geschikt, de repo staat klaar en de eerste stappen voor de analyse zijn duidelijk.

09/09

Plan:

  • Artikel (Parreno et al. 2024) lezen.
  • De belangrijkste onderdelen van het artikel uitwerken.
  • De eerste projecttaken structureren in Trello.

Gedaan:

  • Abstract (Parreno et al. 2024) gelezen.
  • Supplementary Info (inclusief code en data availability) gechecked en gelezen.
  • Trello bullets toegevoegd voor 1. Artikel lezen en 2. Stappen experiment uitlijnen.
  • Beter beeld gekregen van de opzet van het experiment en de data.

Conclusie: Na het lezen en opdelen van de paper in kopjes, is duidelijk geworden wat er voor de presentatie genoemd moet worden. De door de paper gebruikte pipeline is ook duidelijk.

10/09

Presentatie PVA is volgende week dinsdag 10:30.

Plan:

  • Volledig artikel lezen.
  • Ruwe presentatie uitlijnen.
  • Onderzoeksvraag formuleren.

Gedaan:

  • (Globaal) volledig artikel gelezen.
  • Ruwe presentatie uitgelijnd.
    • Dia per rubric blok (“Experiment”, “Experimentele Opzet”, “Materiaal en Methode / Analyse Pipeline”)
    • Criteria in notes v.d. pp. voor gemakkelijke toegang.
    • Steekwoorden en steekzinnen naar wat ik informatief vind:
      • dia 1: Onderzoeksdoel. Context in biologie. Resultaten / conclusies koppelen. What is indel? (insertion/ deletion). What is snv (single nucleotide variant). Miss noemen: grote groepen in de figuur: Downstream, intergenic, intronic -> wat zijn dit allemaal.
      • dia2: Wat voor vliegjes, geslacht. Waar dna vandaan gehaald (eye imaginal disc, wat is dit?). Groepen/ Condities. “thermosensitive ph-RNAi system enabling the reversible KD of ph.”
      • dia3: Referentie genoom: dm6 version of the Drosophila genome.

Conclusie: De presentatie krijgt een eerste opzet. De belangrijkste onderwerpen uit het artikel staan erin, de kerninhoud, begrippen en onderdelen moeten nog verder worden uitgezocht.

11/09

Plan:

  • Rollen verdelen voor PvA pres.

Gedaan:

  • Rollen verdeeld
  • Gestuurd naar team:
Berichten

Dennis: “Hey, zullen we even bespreken wie wat precies gaat doen dit weekend voor de presentatie.”

Jesse: “Yep, is het een goed idee als we allemaal even 1 van de trello PVA punten pakken? Dus als zelf kiezen i guess, voeg dan ff jezelf toe aan de kaart die jij hebt gekozen. Ik doe wel degene die overblijft.”

Jesse: “En dan misschien maandag op school even een onderzoeksvraag voor ons opstellen?”

\(\\\) De genoemde Trello kopjes had ik als volgt aangemaakt:

  • PVA: Pipeline (die van ons en die van hun)
  • PVA: Experiment introductie
  • PVA: Experimentele opzet
  • PVA: Materiaal en methoden

Conclusie: De taken zijn nu verdeeld en de Trello-koppen geven een duidelijke indeling voor de presentatie.

13/09

Plan:

  • PVA: Pipeline is als laatste punt over dus die zal ik uitwerken.
  • Gebruikte pipeline van het onderzoek begrijpen.
  • Alternatieven voor de GATK stappen uitzoeken.
  • Eigen pipeline versie met de GATK stappen geswapped met alternatieve tools.

Gedaan:

  • Pipeline van het onderzoek in tabelvorm in de presentatie gezet.

  • Uitgezocht waar GATK voor werd gebruikt

  • GATK alternatieven gevonden:

    • VarScan2 (Koboldt et al. 2012):
      • Voordeel: Simpel uit te voeren d.m.v .jar bestand.
      • Nadeel: Minder gericht op somatische variatiedetectie en comprimeert de volgorde van de stappen wanneer deze tool gebruikt zou moeten worden.
    • Strelka2 (Kim et al. 2018):
      • Voordeel: Specifiek gericht op somatische variatie detectie, net zoals ook de MuTect2 package voor GATK deed in het originele onderzoek.
      • Nadeel: Meer configuratie en langzamer. \(\\\)

  • Keuze voor onze vervangende pipeline is Strelka2 geworden, ik wil graag de volgorde van stappen zo dicht mogelijk bij het origineel houden en daarin is Strelka2 de beste drop-in replacement (vergelijkend met MuTect2 (Cai et al. 2016) ). Daarnaast zal de snelheid van het programma niet heel veel uitmaken aangezien we op het bin netwerk de analyse mogen doen. :)

  • Tabel van onze versie v.d. pipeline uitgewerkt.

  • Overwegingen en keuzes in de presentatie gezet, volgens de rubric (behalve expliciet referentie genoom informatie, past niet echt meer op de slide dus die moet maar bij materiaal en methode. Staat wel in de tabel maar niet als hoofdpunt zoals in de rubric genoemd).

  • .bib bestand om mijn referenties bij te houden voor dit logboek aangemaakt.

  • Placeholder op de slide voor een terugkoppeling naar onze nog te vormen onderzoeksvraag.

  • Keuze gemaakt om UCSC gene annotation buiten de presentatie te laten, aangezien het een vrij triviale stap is vergeleken met de andere stappen/tools.

Conclusie: De aangepaste pipeline gebruikt Strelka2 als vervanging voor MuTect2. Deze keuze past het beste bij het doel om de volgorde van de oorspronkelijke pipeline zoveel mogelijk te behouden.

14/09

Plan:

  • Data downloaden.
  • Helpen bij de presentatieonderdelen van de anderen, omdat de rest nog niet klaar is.

Gedaan:

  • Accessionlijst gedownload.
  • Het data-downloadbestand klaargemaakt met de code die WERD in Teams heeft gezet.
  • Het downloaden van de data gestart.
  • Een steekwoordenbestand aangemaakt in Teams.
  • De steekwoorden voor mijn dia’s in het bestand gezet.
Steekwoorden slide 6/7

Het originele onderzoek maakte gebruik van de volgende bio-informatische tools: CHECK BORD.

Ze gebruikten hierbij het Drosophila dm6-genoom als referentiegenoom.

SV’s zijn grote DNA-veranderingen, zoals deletions, duplications, insertions, inversions en translocations, meestal groter dan 50 baseparen.

Copy-number variation is een verandering in het aantal gelijke DNA-segmenten in het genoom.

Wij gebruiken geen GATK en ook geen MuTect2. Daarom worden de stappen SNP/InDel calling, variant filtering en somatic variant selection vervangen.

Na het uitvoeren van de pipeline wordt de data in R geanalyseerd.

Onze pipeline verschilt op drie stappen van de originele pipeline. Wij gebruiken geen GATK, maar Strelka2. Strelka2 is een somatische variantcaller die ongeveer vergelijkbaar is met MuTect2 op het gebied van gevoeligheid en nauwkeurigheid. CHECK BORD VOOR REST.

\(\\\)

  • Samen met de groep een voorlopige onderzoeksvraag bedacht: “Kunnen we met de gegeven data bevestigen dat het onderzoek geen carcinogene mutaties bevat?” De vraag is nog niet definitief; de bedoeling is om een simpele en toetsbare vraag voor het project te formuleren. - De data-download is klaar. - Het fasterq-dump-blok gestart. - De download verliep niet helemaal soepel: het script gaf een ncbi_error_file terug en er bleek een probleem te zijn met de tijdelijke opslag. Het fasterq-dump command schreef naar mijn eigen user storage en kon daardoor geen nieuwe dumps meer aanmaken toen deze vol raakte. Dit is opgelost door een nieuwe temp directory aan te maken in dezelfde projectmap op assemblix (/students/2026-2027/Thema05/GenTeam1/temp), zodat de dump genoeg ruimte heeft om te draaien. Daarna is het script aangepast met het nieuwe argument.

09/10: De onderzoeksvraag is aangepast na verduidelijkeing dat met “carcinogene” vaak alleen somatische mutaties bij mensen beschrijft.

Conclusie: De data-download is voorbereid en de voorlopige onderzoeksvraag staat op papier. De tijdelijke opslag was een probleem, maar is opgelost.

15/09

Plan:

  • Een voldoende halen voor de presentatie.

Gedaan:

  • Onze presentatoren extra uitleg gegeven.

  • De GATK-stappen op slide 6 met een rode rand gemarkeerd.

  • Gecontroleerd hoe mijn dia’s aansluiten op de rubric: Referentiegenoom (inclusief versie), analysetools en workflow worden nauwkeurig beschreven. Keuzes in de analyse worden toegelicht en gekoppeld aan de onderzoeksvraag. Begrip van de achterliggende methodologie en mogelijke alternatieven worden getoond.

  • Het referentiegenoom staat bij de mappingstap vermeld als het Drosophila dm6-genoom.

    • Alle analysetools en de volgorde van de stappen worden redelijk nauwkeurig besproken. Ik leg uit wat elke stap doet, maar niet heel gedetailleerd, omdat dit anders niet op de dia past.
    • De keuzes in de analyse worden toegelicht. Ik onderbouw de keuze voor Strelka2 als vervanging voor GATK met informatie uit een onderzoek.
    • Strelka2 wordt gekoppeld aan de onderzoeksvraag. Omdat we de oorspronkelijke workflow zo veel mogelijk willen behouden, is Strelka2 een geschikte vervanging voor MuTect2.
    • VarScan2 wordt als alternatief getoond en er wordt uitgelegd waarom deze tool minder goed bij onze workflow past. \(\\\)

  • Een at-taak ingepland om de fasterq-dump-stap op 16/09 om 01:00 uit te voeren.

Conclusie: Mijn slide sluit op hoofdlijnen aan op de rubric. De workflow, de gemaakte toolkeuzes en het referentiegenoom zijn opgenomen, maar de uitleg op de dia’s blijft bewust kort. De presentatie zelf ging ok.

16/09

Plan:

  • De uitvoer van de at-task controleren.

Gedaan:

  • De uitvoer gecontroleerd; ik kreeg ook al output via e-mail.
  • De documentatie van seqkit gelezen voor het subsetten van fastq.
  • Met seqkit stats statistieken van de gedownloade bestanden opgevraagd. (afgebroken omdat het lang duurde)
  • De versie op de bin-systemen gecontroleerd. Het sample2-commando is daar niet beschikbaar dus een oudere versie.

Conclusie: De taak levert output op en een aantal bestanden zijn met seqkit stats gecontroleerd. Het sample2-commando is niet beschikbaar op de bin-systemen.

21/09

Plan:

  • Een subset maken.

Gedaan:

  • Een subset gemaakt met seqtk (niet met seqkit zoals ik eerder dacht, maar seqtk zoals aanbevolen):
seqtk sample -s 8 fastq/SRR22984519_1.fastq 100000 \
  > subset/SRR22984519_1_subsample.fastq
seqtk sample -s 8 fastq/SRR22984519_2.fastq 100000 \
  > subset/SRR22984519_2_subsample.fastq

  • Gecontroleerd of de read-ID’s in de twee paired-endbestanden overeenkomen. Dat is het geval, dus de subset is goed aangemaakt.

  • Wat permissions in de GenTeam1 folder aangepast:

chmod -R 775 <folders>

\(\\\)

  • fastp uitgevoerd als eerste stap van de pipeline, met de argumenten uit de paper:

  • -g -q 5 -u 50 -n 15 -l 150 –min_trim_length 10 –overlap_diff_limit 1–overlap_diff_percent_limit 10

    • g : Poly-g trimming
    • q : the quality value that a base is qualified
    • u : how many percents of bases are allowed to be unqualified (0~100)
    • n : if one read’s number of N base is >n_base_limit, then this read/pair is discarded.
    • l : reads shorter then lenght will be discarded.
    • min_trim_length : ?
    • overlap_diff_limit : maximum number of mismatched bases to detect overlapped region of PE reads.
    • overlap_diff_percent_limit : maximum percentage of mismatched bases to detect overlapped region of PE reads.

  • Het uiteindelijke command is:

fastp -i SRR22984519_1_subsample.fastq \
  -I SRR22984519_2_subsample.fastq \
  -o SRR22984519_1_clean.fq \
  -O SRR22984519_2_clean.fq \
  -g -q 5 -u 50 -n 15 -l 150 \
  --overlap_diff_limit 1 \
  --overlap_diff_percent_limit 10

  • --min_trim_length is niet meer beschikbaar in versie 0.24.0 en wordt daarom niet gebruikt.

    • Dit heeft gevolgen voor de reproduceerbaarheid, al was reproduceerbaarheid van de paper al een kwestie.

  • Het aantal reads verschilt voor en na QC met ongeveer 25.000 reads.

  • BWA (Burrows-Wheeler Aligner) uitgevoerd met het Drosophila dm6-genoom:

    • “Then, sequencing reads were aligned to the dm6 version of the Drosophila genome using Burrows–Wheeler aligner with default parameters”

    • Het Drosophila dm6-genoom gedownload via UCSC:

wget https://hgdownload.soe.ucsc.edu/goldenPath/dm6/bigZips/dm6.fa.gz
  • BWA op dm6 geindexeerd.
bwa index refgenome/dm6.fa
  • De reads met BWA op dm6 gemapt:
bwa mem refgenome/dm6.fa \
  SRR..._1_clean.fq SRR..._2_clean.fq > SRR...sam

BWA version: 0.7.18-r1243-dirty \(\\\) - De gebruikte commando’s zoals hierboven in een gezamenlijk bestand in Git gezet (used_commands.txt).

  • Het SAM-bestand naar een BAM-bestand geconverteerd en gesorteerd:
samtools sort -o sorted.bam SRR...sam

MarkDuplicates:

  • De Picard-jar gedownload via SourceForge.
  • Geprobeerd Picard MarkDuplicates op het .bam-bestand uit te voeren:
java -jar ~/Downloads/picard.jar MarkDuplicates \
  I=sorted.bam O=marked_dupes.bam M=metrics.txt
  • Dit lukt niet; er verscheen een foutmelding “line doesnt start with numeric”.
    • Uit het debuggen van het SAM- en geconverteerde BAM-bestand met head en less bleek dat het probleem niet per se aan het bestand zelf lag.
    • Een mogelijke oplossing gevonden op het GATK-forum. AddOrReplaceReadGroups met Samtools of Picard lijkt hieruit als oplossing mogelijk te zijn.
  • De fastp-documentatie verder gelezen om uit te zoeken hoe alle 101 FASTQ-bestanden kunnen worden gecontroleerd en getrimd, met 1 HTML-document als uitvoer. Zie Batch processing fastq with fastp.

Conclusie: De subset staat klaar en de eerste stappen van de pipeline zijn uitgevoerd. De read-ID’s komen overeen en het aantal reads daalt met ongeveer 25.000 na fastp. De reads zijn gemapped op dm6 en geconverteerd naar een BAM-bestand. Picard geeft een foutmelding bij het verwerken van het BAM-bestand.

23/09

Plan:

  • Picard debuggen met de eerder gevonden commands.
  • Picard runnen en daarna verdergaan met Strelka2.
  • Proberen fastp op de volledige dataset te runnen.

Gedaan:

  • Het BAM-bestand gecontroleerd op read groups:
samtools view -H SRR_sorted.bam | grep '^@RG'
java -jar ~/Downloads/picard.jar ValidateSamFile \
  I=SRR_sorted.bam MODE=SUMMARY

Er worden geen read groups gevonden. De uitvoer geeft ERROR:MISSING_READ_GROUP terug. Dit bevestigt dus dat er geen RGs aanwezig zijn.

  • Uitgezocht hoe read groups met Picard kunnen worden toegevoegd met AddOrReplaceReadGroups.

  • Opgezocht wat read groups zijn en welke informatie ze bevatten. (Niet veel wijzer van geworden.)

  • Claude opperde dat de library die te gebruiken zou moeten worden voor het toevoegen van Read Groups per sample wel in de Run selector site te vinden was. Daarvoor heb ik het volgende commando uitgevoerd:

curl -o PRJNA918577_runinfo.tsv \
  'https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA918577&result=read_run&fields=study_accession,sample_accession,experiment_accession,run_accession,library_name,library_strategy,library_layout,instrument_model,fastq_ftp&format=tsv'

De output is daarna teruggegeven aan Claude. Hierna werd duidelijk dat de exacte library die in het originele experiment is gebruikt waarschijnlijk niet uitmaakt.

  • Hulp gevraagd om te zien wat we kunnen doen met de AddOrReplaceReadGroups.
  • PicardCommandLine AddOrReplaceReadGroups succesvol uitgevoerd met de volgende instellingen (usage example van AddOrReplaceReadGroups overgenomen):
PicardCommandLine AddOrReplaceReadGroups \
  I=SRR_sorted.bam O=RG_SORTED.bam \
  RGID=4 RGLB=lib1 RGPL=ILLUMINA RGPU=unit1 RGSM=20

\(\\\) - Het BAM-bestand opnieuw gecontroleerd:

PicardCommandLine ValidateSamFile I=RG_SORTED.bam MODE=SUMMARY

De uitvoer geeft aan dat er geen fouten zijn gevonden. - Hierna is MarkDuplicates uitgevoerd:

PicardCommandLine MarkDuplicates \
  I=RG_SORTED.bam O=RG_SORTED_marked_dupes.bam M=metrics.txt
  • fastp mogelijk batch processing commando:
python parallel.py \
  -i ./fastq \
  -o ./clean \
  -r . \
  -p {threads} \
  -a \
  '-g -q 5 -u 50 -n 15 -l 150 --overlap_diff_limit 1 --overlap_diff_percent_limit 10'

Dit gebruikt dezelfde opties als het eerdere commando per sample, maar is bedoeld voor het parallel.py-script van fastp.

  • fastp parallel uitgevoerd met 10 threads.
  • Met -h gecontroleerd hoe het onbekende script paired reads groepeert:
    • parallel.py zoekt standaard naar “R1” en “R2” in de bestandsnamen. Onze samplenamen hebben “_1” en “_2”, waardoor de paired reads niet goed aan elkaar werden gekoppeld.
    • De run opnieuw uitgevoerd met aangepaste argumenten. Omdat de vorige run ongeveer 80 minuten duurde, plan ik deze run weer in met at en het volgende .sh script:
python parallel.py \
  -i ./fastq \
  -o ./clean \
  -r ./reports \
  -p {threads} \
  -1 _1 -2 _2 \
  -a '-g -q 5 -u 50 -n 15 -l 150 --overlap_diff_limit 1 --overlap_diff_percent_limit 10'

Conclusie: Er zijn read groups toegevoegd met AddOrReplaceReadGroups en de BAM-controle geeft geen fouten. Bij de batchverwerking van fastp blijkt de standaardkoppeling van paired reads niet bij onze bestandsnamen te passen, waardoor een nieuwe run nodig was.

26/09

Plan:

  • Het verschil in reads voor en na fastp uitzoeken:
    • De tabel uit reports/overall.html uitlezen.
  • Strelka2 op (mijn i guess?) subset uitvoeren.

Gedaan:

  • De output van de at-run gelezen. De taak gaf python not found terug, omdat ik python i.p.v. python3 in het script had geschreven. Ik heb het script handmatig gestart.

  • Ik heb gekeken naar Strelka2, waarbij ik mij heb gerealiseerd dat we waarschijnlijk een grote fout hebben gemaakt.. We willen uiteindelijk vergelijken tussen de normaal groep en een “gemuteerde” groep. Maar we hebben niet geselecteerd op welke fastq bestanden welke groep betreffen, en ook maar 1 enkel sample i.p.v. een control en een tumor sample..

  • Om er achter te komen welke samples wel gebruikt kunnen worden heb ik de metadata van de accessions gedownload van de SRA Run Selector.

  • Om de goede samples te vinden moeten we zoeken naar “not applicable” of een lege kolom in de chip_antibody kolom. Wanneer deze kolom niet leeg is verwijst het naar een antibody dat gesequenced is door ChIP-seq. Dat is niet wat we nodig hebben voor een whole-genome somatische variantanalyse. Voor onze onderzoeksvraag willen we whole genome samples kunnen vergelijken. Daarnaast moeten de verschillende condities ook worden meegenomen.

  • Geprobeerd de metadata te filteren met:

cat SraRunTable.csv | grep -E 'not applicable|'

  • Dit werkte niet goed voor het zoeken naar lege waarden in chip_antibody.

  • Uiteindelijk een door Claude gegenereerd simpel python bestand om de lijst te subsetten naar “not applicable” en lege chip_antibody rijen gerunned.

    • Hieruit kwam een selectie van mogelijke samples. \(\\\)

  • De accessionnummers in de Run Selector bekeken. Hieruit bleek dat de oorspronkelijke subset één nummer buiten een juiste groep lag. Deze sample kwam uit een specifieke regio en was dus geen whole-genome sample. Vanaf SRR22984520 stonden meerdere whole-genome samples: 20-23 waren Transient RNAi, 24-25 waren Constant RNAi en 26-27 waren Control. Deze samples lijken tot dezelfde batch te behoren (op basis van metadata).

  • Ik kies SRR22984524 als nieuwe tumor-/treated sample en SRR22984526 als control sample.

  • Nu moeten alle uitgevoerde stappen opnieuw worden gedaan voor beide samples.

(Ik heb dit aan/met Dennis en Vani verteld/besproken, maar zij waren nog niet zo ver en waren daarom nog niet helemaal overtuigd. Toen zij bij Strelka2 aankwamen, heeft Dennis Marcel om hulp gevraagd. Daarna werd duidelijk dat samples 24 en 26 inderdaad een juiste keuze waren.)

Uitvoeren stappen op nieuwe data: \(\\\) 1. Seqtk:

seqtk sample -s 31 ../../fastqs/SRR22984524_1.fastq 150000 > subset/Constant_1.fastq
seqtk sample -s 31 ../../fastqs/SRR22984524_2.fastq 150000 > subset/Constant_2.fastq
seqtk sample -s 31 ../../fastqs/SRR22984526_1.fastq 150000 > subset/Control_1.fastq
seqtk sample -s 31 ../../fastqs/SRR22984526_2.fastq 150000 > subset/Control_2.fastq

De bestanden zijn snel gecontroleerd op overeenkomende read-ID’s met head. \(\\\) 2. fastp (voor de volledigheid opnieuw uitgevoerd, ook al was dit met het eerder batch script al gedaan)

fastp -i Constant_1.fastq -I Constant_2.fastq \
  -o Constant_1_clean.fq -O Constant_2_clean.fq \
  -g -q 5 -u 50 -n 15 -l 150 \
  --overlap_diff_limit 1 \
  --overlap_diff_percent_limit 10
fastp -i Control_1.fastq -I Control_2.fastq \
  -o Control_1_clean.fq -O Control_2_clean.fq \
  -g -q 5 -u 50 -n 15 -l 150 \
  --overlap_diff_limit 1 \
  --overlap_diff_percent_limit 10

\(\\\) 3. BWA

bwa index refgenome/dm6.fa

bwa mem refgenome/dm6.fa subset/Constant_1_clean.fastq subset/Constant_2_clean.fastq > bwa_res/Constant.sam
bwa mem refgenome/dm6.fa subset/Control_1_clean.fastq subset/Control_2_clean.fastq > bwa_res/Control.sam

samtools sort -o sorted_constant.bam Constant.sam
samtools sort -o sorted_control.bam Control.sam

\(\\\) 4. Picard

PicardCommandLine AddOrReplaceReadGroups \
  I=sorted_constant.bam O=rg_constant.bam \
  RGID=4 RGLB=lib1 RGPL=ILLUMINA RGPU=unit1 RGSM=20
PicardCommandLine AddOrReplaceReadGroups \
  I=sorted_control.bam O=rg_control.bam \
  RGID=4 RGLB=lib1 RGPL=ILLUMINA RGPU=unit1 RGSM=20
# controle op errors
PicardCommandLine ValidateSamFile I=rg_control.bam MODE=SUMMARY
PicardCommandLine ValidateSamFile I=rg_constant.bam MODE=SUMMARY
PicardCommandLine MarkDuplicates \
  I=rg_constant.bam O=constant_mds.bam M=metrics_const.txt
PicardCommandLine MarkDuplicates \
  I=rg_control.bam O=control_mds.bam M=metrics_control.txt

\(\\\) 5. Strelka2

Strelka2 is nieuw terrein.

  • Strelka2 gedownload en geïnstalleerd via releases op Github.
wget https://github.com/Illumina/strelka/releases/download/v2.9.10/strelka-2.9.10.centos6_x86_64.tar.bz2
tar -xf strelka-2.9.10.centos6_x86_64.tar.bz2

  • BAM bestanden geindexeerd zoals aangegeven in de documentatie van Strelka2.

    • dm6 geindexeerd.
    samtools index constant_mds.bam
samtools index control_mds.bam
samtools faidx refgenome/dm6.fa
  • Structural variant caller Manta maakt Strelka preciezer, dus die wil ik toch ook uitvoeren.
  • Voor beide programmas blijkt python2 nodig volgens de documentatie. Deze versie is niet aanwezig op assemblix, dus ik wil proberen via bioconda de beide programmas te downloaden.
    • Dit betekent dat de eerdere gedownloade tar files verwijderd kunnen worden. \(\\\)
  • Bioconda installatie:
conda create -n variantcalling \
  -c conda-forge -c bioconda strelka manta
  • Installeert Bioconda, conda-forge en strelka en manta in een environment genaamd “variantcalling”.

Manta-configuratie:

conda activate variantcalling

configManta.py \
  --normalBam ../picard_res/control_mds.bam \
  --tumorBam ../picard_res/constant_mds.bam \
  --referenceFasta ../refgenome/dm6.fa \
  --runDir ../manta_out
  • Dit geeft manta_out/runWorkflow.py terug, waarmee de workflow kan worden gestart.
  • runWorkflow.py heeft een -j optie voor het aantal threads; deze zet ik op 15.
    • De output bevat drie folders (evidence, variants en stats). In variants staan onder andere candidateSV.vcf.gz, candidateSmallIndels.vcf.gz, diploidSV.vcf.gz en somaticSV.vcf.gz.

Strelka2-configuratie:

configureStrelkaSomaticWorkflow.py \
  --normalBam picard_res/control_mds.bam \
  --tumorBam picard_res/constant_mds.bam \
  --referenceFasta refgenome/dm6.fa \
  --indelCandidates manta_out/results/variants/candidateSmallIndels.vcf.gz \
  --runDir strelka_out
  • Dit geeft strelka_out/runWorkflow.py terug, waarmee de workflow kan worden gestart.
  • runWorkflow.py heeft een -j optie voor het aantal threads; deze zet ik op 15.
  • De workflow geeft errors. Zelfs het uitvoeren van --version geeft een segmentation fault. Ik heb Claude gevraagd om opties, waarna het idee om een docker/podman versie van Strelka2 te runnen naar voren kwam.
  podman pull quay.io/biocontainers/strelka:2.9.10--h9ee0642_1
  podman run --rm \
    quay.io/biocontainers/strelka:2.9.10--h9ee0642_1 \
    configureStrelkaSomaticWorkflow.py --version
  • Dit geeft versie 2.9.10 terug. Daarmee lijkt Strelka2 in de container wel te werken.

Hoe nu Strelka te configureren:

podman run --rm \
  -v /students/2026-2027/Thema05/GenTeam1/andere_subsets/j_ss:/data \
  quay.io/biocontainers/strelka:2.9.10--h9ee0642_1 \
  configureStrelkaSomaticWorkflow.py \
    --normalBam /data/picard_res/control_mds.bam \
    --tumorBam /data/picard_res/constant_mds.bam \
    --referenceFasta /data/refgenome/dm6.fa \
    --indelCandidates /data/manta_out/results/variants/candidateSmallIndels.vcf.gz \
    --runDir /data/strelka_out
  • --rm optie: verwijderd de container meteen na het runnen weer.
  • -v optie: volume waarin de container zal gaan runnen; het pad voor ‘:’ is de standaard folder die in de container gekopieerd zal worden in de folder na ‘:’. \(\\\) En de workflow:
# -m local is voor lokaal runnen, strelka heeft de optie voor remote
podman run --rm \
  -v /students/2026-2027/Thema05/GenTeam1/andere_subsets/j_ss:/data \
  quay.io/biocontainers/strelka:2.9.10--h9ee0642_1 \
  /data/strelka_out/runWorkflow.py -m local -j 15

(Podman en samtools werken overigens niet op de assemblix-server. Deze commandos zijn via SSH op NUC 412 uitgevoerd.)

Beide commando’s lijken te werken en er staan resultaten in de outputmap. Voor de verdere analyse is het belangrijk om te controleren welke varianten PASS in de FILTER-kolom hebben. In de paper staat hierover: “Only SNVs and InDels variants that passed Mutect2 filtering (FILTER=”PASS”) were considered for downstream analyses.”

Conclusie: De samplekeuze staat nu op SRR22984524 als tumor en SRR22984526 als control. De eerdere stappen worden voor deze samples opnieuw uitgevoerd. Lokaal geeft Strelka2 nog een segmentation fault, terwijl de configuratie en workflow in de container vooralsnog werken. Samtools en Podman werken niet op assemblix, wel op de NUCs.

28/09

Plan:

  • CNVnator en Breakdancer uitvoeren.

Gedaan:

conda create -n cnvnator -c conda-forge -c bioconda cnvnator

\(\\\) - De coverage gecontroleerd met Samtools. CNVnator heeft voldoende coverage nodig voor een bruikbaar resultaat.

samtools coverage tumor.bam

~0.3x

\(\\\) - Breakdancer is ook via Bioconda te installeren:

conda create -n bdancer -c conda-forge -c bioconda breakdancer

\(\\\) - Het verschil in reads tussen voor en na QC van de subsets gecontroleerd:

  • Constant: voor fastp 300.000 reads en 45 miljoen basen; na fastp 295.450 reads en 44,3 miljoen basen.

  • Control: voor fastp 300.000 reads en 45 miljoen basen; na fastp 296.404 reads en 44,4 miljoen basen. \(\\\)

  • CNVnator heeft meerdere voorbereidende stappen nodig, waaronder het opsplitsen van dm6.fa in individuele chromosomen.

  • Tussendoor aan het team doorgegeven dat ik samples SRR22984524 en SRR22984526 al had gekozen.

    • Ook aangegeven dat ik Manta voor Strelka2 had uitgevoerd. \(\\\)

  • Een Galaxy-workflow gestart met de tools uit de pipeline om de analyse later beter reproduceerbaar te maken.

    • Galaxy kan sommige fastp-opties niet uitvoeren (–overlap_diff_limit, –overlap_diff_percent_limit). Op de subset scheelt dit ongeveer 200.000 basen met het wel complete fastp-command. Misschien is het verschil op het complete sample minder groot. \(\\\)

  • Op dit moment bevat de Galaxy-workflow de stappen van sample inladen tot en met Strelka2. Breakdancer en CNVnator lijken niet als tools beschikbaar te zijn in Galaxy, dus hiervoor moet mogelijk een alternatief worden gezocht.

    Galaxy Workflow tot zo ver
    Galaxy Workflow tot zo ver

Conclusie: De read-aantallen dalen na fastp met ongeveer 4.500 reads voor Constant en 3.600 reads voor Control. De coverage van de subsets lijkt te beperkt voor een betrouwbare CNV-analyse, dus de volledige samples zijn waarschijnlijk geschikter. De Galaxy-workflow staat voorlopig tot Strelka2, sommige fastp opties zijn niet beschikbaar in Galaxy.

30/09

Plan:

  • CNVnator en Breakdancer uitvoeren.
  • Kijken of er resultaten uit de subset komen.
  • Alternatieven voor CNVnator en Breakdancer zoeken voor de Galaxy-workflow.
  • De pipeline op een volledig sample uitvoeren.

Gedaan:

  • CNVnator uitgevoerd op de Control- en Constant-subset via Conda.
  • De read mapping uit het BAM-bestand gehaald:
cnvnator -root control_out.root -tree control_mds.bam

\(\\\) - dm6.fa opgesplitst in individuele chromosomen voor het read-depth histogram:

cd chromosomes
csplit -s -z /students/2026-2027/Thema05/GenTeam1/refgenome/dm6.fa '/>/' '{*}'
for i in xx* ; do \
  n=$(sed 's/>// ; s/ .*// ; 1q' "$i") ; \
  mv "$i" "$n.fa" ; \
 done

Via: housegordon

  • Eerste CNVnator-run met een bin size van 500 uitgevoerd:
cnvnator -root control_out.root -his 500 -d chromosomes/
cnvnator -root control_out.root -stat 500
  • De gerapporteerde read depth is laag (\(xe^{-13}\)). Met samtools stats komt de coverage uit op ongeveer 0.3x. Daarom opnieuw gerekend met alleen de belangrijkste chromosomen aangezien de random/un chromosomen weinig data hadden:
cnvnator -root control_out.root -chrom chr2L chr2R chr3L chr3R chr4 chrX -tree control_mds.bam
cnvnator -root control_out.root -chrom chr2L chr2R chr3L chr3R chr4 chrX -his 500 -d chromosomes/
cnvnator -root control_out.root -chrom chr2L chr2R chr3L chr3R chr4 chrX -stat 500

  • De read depth blijft laag. Met een bin size van 10.000 komt de read depth rond 20 uit.

  • Met een bin size van 18.000 komen de ratios volgens cnvnator -eval uit op ongeveer 3,91 en 4,20 voor Constant en 4,01 en 4,12 voor Control. Dit ligt rond de ratio van ongeveer 4 die in de CNVnator-documentatie/Github Issues wordt genoemd en aanbevolen.

  • Nu kunnen de laatste stappen worden uitgevoerd:

cnvnator -root control_out.root -chrom chr2L chr2R chr3L chr3R chr4 chrX -partition 18000
cnvnator -root control_out.root -chrom chr2L chr2R chr3L chr3R chr4 chrX -call 18000 > control_cnv.txt
  • De output van Constant is leeg, en de output van Control bevat slechts 1 CNV. Omdat de subset maar beperkte coverage heeft, is het waarschijnlijk betrouwbaarder/meer waardevol om op de gehele samples te testen.

Control_cnv.txt:

CNV_type coordinates CNV_size normalized_RD e-val1 e-val2 e-val3 e-val4 q0
deletion chr2L:23130001-23346000 216000 0.631105 0.00969536 3396.56 0.00969536 3396.56 0.437063

Breakdancer

  • Een configuratiebestand gegenereerd met bam2cfg.pl volgens de documentatie op Github:
# zelfde voor constant
bam2cfg.pl -g -h picard_res/control_mds.bam > control.cfg

“If you have a single bam file that contains multiple libraries, make sure that the readgroup and library information are properly encoded in the sam/bam header, and in each alignment record, otherwise bam2cfg.pl may fail to produce a correct configuration file.”

Dit heeft misschien ook invloed aangezien ik de readgroup en library informatie handmatig heb toegevoegd.

  • control.cfg is leeg. Daarom alleen geprobeerd om met Constant verder te gaan.
breakdancer-max constant.cfg > constant_sv.txt

Hier kwam de volgende informatie uit:

Chr1 Pos1 Orientation1 Chr2 Pos2 Orientation2 Type Size Score Num Reads Num Reads (lib) Sample
chr3L 11611403 2+0- chr3L 11611992 0+2- DEL 450 99 2 2 constant_mds.bam
chr3L 17799698 2+0- chr3L 17801596 0+3- DEL 1746 99 2 2 constant_mds.bam
chrM 1 0+15- chrM 19449 13+0- ITX 18878 99 12 12 constant_mds.bam
chrUn_DS485995v1 1 50+23- chrUn_DS485995v1 160 50+23- ITX -248 99 21 21 constant_mds.bam
chrUn_DS484226v1 1 85+66- chrUn_DS484226v1 32 85+66- ITX -224 99 18 18 constant_mds.bam
chrUn_DS484226v1 27 85+66- chrX 22433013 25+27- CTX -327 99 7 7 constant_mds.bam
chrUn_DS485255v1 563 11+11- chrUn_DS485255v1 678 11+11- ITX -215 99 8 8 constant_mds.bam
chrUn_DS485249v1 144 2+2- chrUn_DS485249v1 205 2+2- ITX -193 99 2 2 constant_mds.bam

Statistics:

Metric Value
Mean insert size 328.48
Standard deviation 78.88
Upper cutoff 730.92
Lower cutoff 113.79
Read length 149.63
Library lib1
Reference length 139788235
Sequence coverage 0.264268
Physical coverage 0.290071

En de histogrammen (control, constant):

Control histogram Constant histogram

Met geen goede control om tegen te vergelijken is er niet veel uit deze resultaten te halen. Een analyse met volledige samples lijkt nu een betere optie. Voor de Galaxy-workflow ga ik alternatieven voor CNVnator en Breakdancer zoeken:

Origineel Alternatief
CNVnator CNVkit
Breakdancer Delly

  • In Galaxy is de CNVkit batch-node gebruikt. De Control- en tumordelen van de workflow zijn met gekleurde vakken gemarkeerd voor duidelijkheid. De bin size voor CNVkit kan via een aparte input-node worden aangepast.

  • CNVkit vraagt ook bij whole-genome data om een BED-bestand. Dit bestand lijkt niet direct in Galaxy gegenereerd te kunnen worden, dus het moet apart worden gemaakt en geüpload. Hier is ook een aparte input-node voor gemaakt.

  • De Galaxy-workflow wordt steeds groter, maar is overzichtelijker dan de tools handmatig uitvoeren:

    Workflow nu
    Workflow nu

De volledige, opgeschoonde (met fastp) bestanden van samples 24 en 26 uit de clean/-map gekopieerd:

cp ../../clean/SRR22984524* .
cp ../../clean/SRR22984526* .

  • De samples handmatig geüpload, omdat de SRA-downloadnode in Galaxy een collectie teruggeeft en het niet heel makkelijk is hoe de juiste samples aan de tumor- en controlkant van de workflow moeten worden gekoppeld.

  • De workflow wordt als laatste met deze samples uitgevoerd.

  • De workflow is gecrashed omdat de aangeleverde samples al door de QC heen waren geweest en dus niet wat fastp verwachtte. Ik ga originele samples uploaden.

Conclusie: De subset geeft ongeveer 0,3× coverage. CNVnator geeft daardoor waarschijnlijk geen betrouwbare uitkomst: Constant levert geen CNV op en Control één. Een analyse met volledige samples lijkt daarom geschikter, hetzelfde geld voor Breakdancer. Voor Galaxy zijn CNVkit en Delly als alternatieven gekozen, en zijn de volledige originele samples nodig.

02/10

Plan:

  • De workflow met de volledige samples uitvoeren.

Gedaan:

  • De originele samples uit de FASTQ-map op /students/2026-2027/Thema05/GenTeam1/fastqs gekopieerd, individueel met gzip gecomprimeerd en als collectie in een tar-bestand gezet voor de upload.
gzip *.fastq
tar -cf samps.tar SRR22984524_1.fastq.gz SRR22984524_2.fastq.gz SRR22984526_1.fastq.gz SRR22984526_2.fastq.gz

  • De naar Galaxy geuploade tar geunzipped met een ingebouwde tool.

  • De workflow uitgevoerd met Control R1: *26_1, Control R2: *26_2, enzovoort.

  • De workflow is gecrasht; CNVkit en Delly geven allebei fouten.

  • De stderr van CNVkit geeft het volgende aan:

Summary: #bins=7046, #reads=0, mean=0.0000, min=nan, max=nan
Skip processing normal.bam with empty regions file ./capture.antitarget.bed

  • Deze melding verschijnt ook voor het tumor-BAM-bestand.

  • De ingebouwde Galaxy tool samtools idxstats uitgevoerd op de input van CNVkit en Delly, dus op de output van MarkDuplicates:

    • Dit geeft ook 0 reads aan. Daarom ook de output van BWA met samtools gecontroleerd om te kijken of daar wel reads in staan.
    • Ook BWA-MEM bevat geen reads. De BWA-input bestaat uit dataset 108, de dm6-index en de output van fastp. De fastp-output lijkt dus niet goed te zijn:

fastp output \(\\\)

  • De input-fastq.gz-collectie lijkt wel goed te zijn, maar bij het inlezen door fastp gaat het mis.
  • De fastp-output geeft aan dat het FASTQ-bestand mogelijk ongeldig is en adviseert om de laatste regels met tail te controleren. Daarom de laatste 20 regels van een sample bekeken met de Galaxy-tool Select Last (lines from a dataset). Deze regels lijken normaal:

Laatste regels \(\\\) - fastp verwacht een fastqsanger-datatype, terwijl de bestanden als fastq.gz zijn aangeleverd. - Het datatype in Galaxy aangepast via Edit Attributes, van fastq.gz naar fastqsanger.gz.

  • fastp los van de Workflow uitgevoerd op 2 samples. De output lijkt nu weer te kloppen:

fastp output los van workflow \(\\\) De volledige workflow wordt opnieuw uitgevoerd om te controleren of deze het nu wel doet.

Conclusie: De eerste workflowrun werkt nog niet, omdat de CNVkit- en Delly-stappen geen reads ontvangen. De oorzaak lijkt bij de fastp-input of het datatype te liggen. Na het aanpassen naar fastqsanger.gz lijkt de fastp-output weer goed, maar de volledige workflow moet nog opnieuw worden uitgevoerd.

04/10

Plan:

  • De logboekstructuur verbeteren en de inleiding duidelijker maken.
  • De spelling, grammatica en formulering van de eerdere dagen verbeteren.
  • Resultaten op een overzichtelijke manier presenteren in tabellen in het logboek.
  • De wederom gecrashte Galaxy Workflow debuggen.

Gedaan:

  • De logboekstructuur verbeterd en de inleiding duidelijker gemaakt.

  • De spelling, grammatica en formulering van de eerdere dagen verbeterd.

  • Resultaten op een overzichtelijke manier gepresenteerd in tabellen in het logboek.

  • De Galaxy-workflow gaf nog steeds fouten in Delly en Strelka2:

    • Delly: [E::bcf_hdr_add_sample_len] Duplicated sample name '20'
    • Strelka2: [E::get_intv] Failed to parse TBX_VCF \(\\\)

  • De Delly-fout leek voort te komen uit duplicate read-group namen. Delly gebruikt de input van beide MarkDuplicates-stappen en daarmee ook beide AddOrReplaceReadGroups-stappen, waarin de BAM-bestanden dezelfde argumenten voor het toevoegen van de Read Groups kregen. Ik heb daarom de read-group-instellingen voor beide samples aangepast naar duidelijke, unieke waarden:

RG_ID RG_SM LB PL PU Opmerking
4 20 lib1 ILLUMINA unit1 Default settings
12 tum tumor1 ILLUMINA unitum Tumor argumenten
6 ctrl control1 ILLUMINA unitctr Control argumenten

  • De Strelka2-fout is waarschijnlijk veroorzaakt door Manta. Strelka2 gebruikt het argument --indelCandidates, dat een outputbestand van Manta verwacht. In de workflow had ik het verkeerde bestand (somaticSV) gekoppeld in plaats van het bestand candidateSmallIndels.vcf.gz. Met deze verandering zou de stap Manta -> Strelka opgelost moeten zijn.

  • De workflow is opnieuw gestart met de aangepaste instellingen.

  • Logboek source vergeleken met de knit.

    • Output knit HTML thema veranderd naar cosmo met bootswatch. M.b.v. bslib::bs_theme_preview().
    • Dark mode knop toegevoegd met bslib: bslib::input_dark_mode()
    • Extra ruimte toegevoegd waar de knit ruimte weghaalde met LaTeX: $\\$

Conclusie: De problemen van deze dag lagen in de read-group-naamgeving en de verkeerde koppeling van Manta-output naar Strelka2. Nadat deze problemen zijn opgelost is de workflow opnieuw gestart. De structuur van het logboek is daarnaast ook veel duidelijker geworden.

05/10

Plan:

  • De Manta -> Strelka2 koppeling gaat nog niet goed in de workflow, dit moet gefixed worden.

Gedaan:

  • Manta’s output is vcf_bgzip. Het datatype veranderen in Galaxy verandert niks aan het bestand zelf:
file *.vcf
# geeft bgzip terug
  • Uncompress node werkt ook niet om de vcf netjes uit de bgzip te halen.
  • Unzip node werkt ook niet.
  • bcftools view werkt wel! Met het output type gezet als uncompressed vcf is de file leesbaar voor Strelka.
CHROM POS ID REF ALT QUAL FILTER
chr2L 5372 . T A . LowEVS
chr2L 5390 . T A . LowEVS
chr2L 5598 . C G . LowEVS
chr2L 5762 . T C . LowEVS
chr2L 5904 . C A . LowEVS
chr2L 5974 . C T . LowEVS
chr2L 7039 . A T . LowEVS
chr2L 7088 . A T . LowEVS
chr2L 7556 . T C . PASS
chr2L 7902 . A G . LowEVS
chr2L 8263 . G A . LowEVS
chr2L 12275 . G C . LowEVS
chr2L 12460 . T C . PASS

Strelka2 SNV resultaat. (Klein overzicht, zie ook de “PASS” waarden welke nodig zijn voor de analyse.) \(\\\)

CHROM POS ID REF ALT QUAL FILTER
chr2L 10588 . AT A . LowEVS
chr2L 18734 . TATA T . LowEVS
chr2L 39030 . T TAAGAG . PASS
chr2L 39031 . CT C . PASS
chr2L 39036 . A AGAGCGG . PASS
chr2L 39958 . T TTCACG . LowEVS
chr2L 39964 . GATGTGGAAAAA G . LowEVS
chr2L 39977 . ACAAG A . LowEVS
chr2L 47084 . TTCTA T . LowEVS
chr2L 54453 . TATAATATATAATA T . PASS
chr2L 55718 . TA T . LowEVS
chr2L 57619 . A ACT . LowEVS
chr2L 60531 . ATTT A . LowEVS

Strelka2 Indels resultaat. (Klein overzicht, zie ook de “PASS” waarden welke nodig zijn voor de analyse.)

  • Enige punt wat ik nog ben vergeten: Ik heb bij CNVkit de bin_size op 19000 gehouden al is de readdepth van de gehele samples natuurlijk groter dan bij de subsets. Ik verwijder de input aan CNVkit met de optie om zelf een bin_size op te geven, CNVkit kan dit zelf uitrekenen op basis van de readdepth.

  • De gemiddelde readdepth van de gehele sample (getest op 1 bam, alleen nucleaire chromosomen) is ongeveer 109x:

samtools coverage *.bam
#rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq
chr2L 1 23513712 17474577 23481327 99.8623 109.835 35.6 55.3
chr2R 1 25286936 18687851 25237357 99.8039 109.152 35.6 52.6
chr3L 1 28110227 20989336 27941302 99.3991 110.237 35.6 53.4
chr3R 1 32079331 23484246 32002751 99.7613 108.337 35.6 55
chr4 1 1348131 1083772 1325139 98.2945 119.041 35.5 53.6
chrX 1 23542271 17823084 23428270 99.5158 111.25 35.5 52.2

Resultaat samtools coverage \(\\\) - Ter vergelijking was de readdepth van de subsets 0.3x..

Conclusie: De workflow loopt nu goed van eind tot eind. De resultaten van de workflow moeten gefilterd en geannoteerd worden.

07/10

Plan:

  • Output vcfs filteren.
  • Begin analyse met R.

Gedaan:

  • CNVkit call node toegevoegd na de cnvkit batch: Deze kiest daadwerkelijk de juiste cnvs.
  • Handmatig Delly output omgezet in vcf met bcftools view.
  • Ik wil een soort template Rmd hebben zodat ook andere samples meteen geanalyseerd kunnen worden.
    • Ik zoek wat meer op over Rmd syntax. Zoals params en andere config opties.
    • Eerste filtering code voor Strelka Delly en CNVkit gemaakt voor filteren op “PASS”.

Conclusie: Eerste stappen voor de analyse gedaan. En (nu echt) allerlaatste probleem in de Workflow verholpen.

08/10

Plan:

  • Analyse afmaken.

Gedaan:

  • Waarden van filteren uit de paper gehaald:
Filterwaarde Source
min_vaf: 0.2 allelic fraction greater than 0.2
sv_min_support: 5 supported by at least five reads
cnv_gain: 0.585 CNVs ..allelic fraction > than 1.5 (\(\log_2(1.5)= 0.585\))
cnv_loss: -0.6 smaller than 0.66 (\(\log_2(0.66)\approx -0.6\))

*“we only retained SNVs or InDels with an allelic fraction greater than 0.2, structural variants that were supported by at least five reads and CNVs with an allelic fraction bigger than 1.5 (duplication) or smaller than 0.66 (deletion).” (Parreno et al. 2024)

Filterwaarden Reden
min_depth: 30 Minimale readdepth voor precizie.
max_normal_alt: 0 0 somatische reads in normaal.
cnv_min_probes: 20 Noise filter voor cnv bins
  • Gebruik variantannotation om vcfs te lezen.
File Source
Strelka_som_snvs.vcf Strelka2 SNV output
Strelka_som_indels.vcf Strelk2 Indel output
manta_score_filtered_variants.vcf Manta variants file (niet indel_candidates voor strelka!)
Delly.vcf Delly SV call output
CNVkit_CALL_manual.cns CNVkit call output (gecentered met -0.082)
CNVkit_segs.cns Ongecallde (is dit een woord: nu wel) segments van CNVs
CNVkit_CALL.cns Oude CNVkit call output

  • Voor de R code heb ik AI als hulpmiddel. Alle code die door AI is geschreven is gelezen/gechecked/aangepast/commented.

  • De resultaten spreken niet voor zich, ik heb er nog tekst bij geschreven.

  • Iframe voor de Galaxy workflow toegevoegd.

    • Dit werkt niet helemaal lekker, ligt waarschijnlijk aan file:// en bioinf.nl/view?. Hopelijk werkt het wel op Github.

  • Onderscheid tussen code-blokjes en outputblokjes duidelijker gemaakt in Dark mode met css styling.

  • VAF resultaten van de paper nog eens opgezocht:

    • “Moreover, 92.8% of the identified SNVs or InDels had an allele frequency below 0.2, precluding them from driving whole-tissue tumours”

  • CNVkit Call opnieuw uitgevoerd met een handmatige median shift van -0.082, de logaritmische waarden in de cns output was niet gecenterd rond 0. Met de handmatige setting zou dit verholpen moeten zijn.

  • code-folding optie toegevoegd aan de YAML-header van dit bestand: maakt het inklappen van code blokjes mogelijk.

  • self_contained optie toegevoegd aan de YAML-header van dit bestand: maakt figuren base64 strings i.p.v. een referentie naar een bestand.

  • Conclusie en beperkingen van de analyse opgeschreven.

  • GenTeam1 folder op BIN netwerk opgeruimd

    • .old folder met eerste (foutieve) gezamenlijke subset resultaten.
    • ge-QCde fastq dir clean/ verwijderd.
    • scripts zoals parallel.py en fadump.sh naar een directory used_scripts.
    • SRA folder verwijderd.

  • j_ss subset folder opgeruimd.

    • contig debugging files verwijderd.
    • Rest zag er wel overzichtelijk uit, heb nog een file toegevoegd die de inhoud van elk directory uitlegd wanneer ls word gerunned.

  • Typefouten uit logboek gehaald.

    • 29 occurrences van FastP of Fastp veranderd naar fastp.

Conclusie: De analyse is klaar. De GenTeam folder is opgeruimd net zoals mijn persoonlijke subset folder. Details zijn nog weggewerkt in dit logboek.

09/10

Plan:

  • n.v.t.

Gedaan:

  • Logboek extra aangepast.
  • Marcel stelde SnpEff voor als resultaat, maar ik heb al een analyse gedaan en ik vond daar juist in dat annoteren van genen niet veel toevoegd voor resultaat en helemaal niks voor het beantwoorden van de onderzoeksvraag. Dus dit ga ik niet meer doen.
  • Code blokjes in analyse voorzien van class.source = "fold-hide".
  • data-external = "1" toegevoegd aan Galaxy embed.
  • offset_path optie toegevoegd aan de config voor gemakkelijker runnen.
  • Fix: Verkeerd percentage en conclusie bij SNVs.

Conclusie: Ik ben (er) klaar (mee). \(\\\) \(\\\)

Analyse

In dit deel van het logboek staat alle code voor het analyseren, filtreren en visualiseren van de output van de pipeline.

Per sectie is uitgelegd wat de code doet, wat de data vertelt en wat voor conclusies uit de resultaten te behalen vallen.

De code blokjes zijn uitklapbaar met de knop.

Om exact dezelfde analyse uit te voeren dien je de code uit dit bestand uit te voeren in GenTeam1/andere_subsets/j_ss/FINAL, pas zo nodig de variabele offset_path aan in de config hieronder.

Config

Dit blok code zorgt ervoor dat belangrijke globale argumenten beschikbaar en gemakkelijk aan te passen zijn.

De bestanden die worden geladen zijn beschikbaar in GenTeam1/andere_subsets/j_ss/FINAL

# Load samples

sample_str <- "24-26" # om mogelijke runs met andere samples mogelijk te maken via een andere map
offset_path <- "../../"
result_dir <- paste(offset_path, "files/", sample_str, sep="")

strelka_snv <- paste(result_dir, "/Strelka_som_snvs.vcf", sep="")
strelka_indel <- paste(result_dir, "/Strelka_som_indels.vcf", sep="")
manta_vcf <- paste(result_dir, "/manta_score_filtered_variants.vcf", sep="")
delly_vcf <- paste(result_dir, "/Delly.vcf", sep="")
cnvkit_cns <- paste(result_dir, "/CNVkit_CALL_manual.cns", sep="") # handmatig geshift bestand

main_chromosomen <- c("chr2L", "chr2R", "chr3L", "chr3R", "chr4", "chrX")

# extra filter argumenten
min_vaf <- 0.2
sv_min_support <- 5
cnv_gain <-      0.585  # 2log(1.5)
cnv_loss <-      -0.6   # ongeveer 2log(0.66)

min_depth <-     30
max_normal_alt <- 0
cnv_min_probes <- 20


library(VariantAnnotation) # vcfs inlezen
library(dplyr)
library(tidyr)
library(ggplot2)

Vcfs inlezen

Wat de code doet: De functie read_strelka() haalt per variant op hoeveel reads de variant dragen in de tumor en in de controle, en rekent daaruit de VAF van de tumor uit. De variant wordt “gehouden” (keep) als hij aan al deze eisen voldoet:

Regel Waarde Reden
filter == "PASS" – De caller zelf keurt de variant goed.
Hoofdchromosomen chr2L, chr2R, chr3L, chr3R, chr4, chrX Losse scaffolds (chrUn) en chrM zijn onbetrouwbaar, chrY is mannelijk en niet van toepassing.
min_depth >= 30 30 reads Een VAF uit een paar reads is onbetrouwbaar.
max_normal_alt <= 0 0 reads De variant mag niet in de controle zitten.
min_vaf >= 0.2 0,2 Drempel uit het artikel (Parreno et al.). Een mutatie die een tumor aanstuurt, moet in 20% v.d. reads zitten. 
# functie die alleen het eerste alt allel pakt: eg. G,T -> G
first_alt <- function(vcf) {
  a <- alt(vcf)
  flat <- as.character(unlist(a))
  starts <- cumsum(c(1, lengths(a)))   # c(1, lengths(a)) gives the offsets;
  first_indices <- starts[seq_along(a)]
  # pull out the first item of each list element
  result <- flat[first_indices]
}

who <- function(vcf) { # identificeerd tumor sample als die begint met 't'
  s <- samples(header(vcf)) # samples returns a character vector of samples (VariantAnnotation)
  t <- grep("^t", s, ignore.case = TRUE) # index van tumor sample rij
  c(tumor = s[t], normal = s[-t]) # indetify the tumour sample, other is normal
}

read_strelka <- function(path, kind) {
  vcf <- readVcf(path, genome = "dm6")
  s <- who(vcf)
  g <- geno(vcf) # genotype information from a VCF file (VariantAnnotation)

  ref <- as.character(rowRanges(vcf)$REF)
  alt <- first_alt(vcf)

  if (kind == "SNV") {  # snvs vcf heeft een andere structuur; snv bewijs per nucleotide
    counts <- function(sample) {
      cbind(g$AU[, sample, 1], g$CU[, sample, 1],
            g$GU[, sample, 1], g$TU[, sample, 1]) # bouw een matrix met row = variants, cols = nucs en x varianties
    }

    i <- seq_along(ref)
    t <- counts(s["tumor"])[cbind(i, match(alt, c("A","C","G","T")))] # tumor alt reads
    r <- counts(s["tumor"])[cbind(i, match(ref, c("A","C","G","T")))] # tumor ref reads
    n <- counts(s["normal"])[cbind(i, match(alt, c("A","C","G","T")))] # normal alt reads
  } else {
    t <- g$TIR[, s["tumor"], 1] # tumor alt/indel reads
    r <- g$TAR[, s["tumor"], 1] # tumor reference reads
    n <- g$TIR[, s["normal"], 1] # normal reads supporting the indel
  }

  tibble(
    kind = kind,
    chrom = as.character(seqnames(vcf)),
    pos = start(vcf),
    filter = as.character(rowRanges(vcf)$FILTER),
    t_depth = t + r,
    t_alt = t,
    t_vaf = t / (t + r),
    n_alt = n
  ) # return data as table
}

# Combine SNV and indel calls into one table
small <- bind_rows(
  read_strelka(strelka_snv, "SNV"),
  read_strelka(strelka_indel, "indel")
)

chrom_len <- seqlengths(
  readVcf(strelka_snv, genome = "dm6")
)[main_chromosomen]

head(small)
## # A tibble: 6 × 8
##   kind  chrom   pos filter t_depth t_alt t_vaf n_alt
##   <chr> <chr> <int> <chr>    <int> <int> <dbl> <int>
## 1 SNV   chr2L  5372 LowEVS     104    16 0.154     5
## 2 SNV   chr2L  5390 LowEVS     108    16 0.148     7
## 3 SNV   chr2L  5598 LowEVS     107    17 0.159    11
## 4 SNV   chr2L  5762 LowEVS     100    20 0.2      12
## 5 SNV   chr2L  5904 LowEVS      80    15 0.188     9
## 6 SNV   chr2L  5974 LowEVS      85    15 0.176     8

Strelka PASS filtering

Strelka2 resultaten filteren met agressievere thresholds.

small <- small %>% mutate(
  pass = filter == "PASS" & chrom %in% main_chromosomen,
  keep = pass & coalesce(t_depth >= min_depth &
                         n_alt <= max_normal_alt & t_vaf >=min_vaf, FALSE))

small %>% group_by(kind) %>%
  summarise(num = n(), PASS = sum(pass), na_overige_filters = sum(keep)) %>%
  dplyr::tibble()
## # A tibble: 2 × 4
##   kind     num  PASS na_overige_filters
##   <chr>  <int> <int>              <int>
## 1 SNV   126297 10025                 65
## 2 indel   9400   593                 12
small %>%
  filter(kind == "SNV", pass) %>%
  summarise(PASS = n(),
            in_nrml = sum(n_alt >= 1),
            percentage = round(100 * mean(n_alt >= 1), 1)) %>%
  knitr::kable()
PASS in_nrml percentage
10025 7989 79.7

Resultaat: Van 126.297 SNV-records houdt Strelka2 er 10.025 over als PASS geclassificeerd. Na de aanvullende filtering blijven er 65 over. Het grootste verlies zit bij de eis dat er geen alt-reads in de controle mogen zitten voor een echte somatische call: 79,7% van de PASS-SNV’s heeft minstens 1 alt-read in de controle. Dat betekent dat de meeste “somatische” PASS-calls waarschijnlijk niet specifiek voor de tumor zijn, al is dit een sterke filter. \(\\\) \(\\\)

Density grafiek: Op de x-as staan de chromosomen, op de y-as het aantal PASS-SNV’s per miljoen baseparen (Mb). Delen door de chromosoomlengte is nodig, anders lijkt een lang chromosoom vanzelf meer te hebben. Elke staaf is gestapeld: blauw is “gehouden”, oranje is “gefilterd”. Omdat er maar 65 gehouden SNVs zijn van de 10.000+ PASS-SNVs, is het blauwe deel nauwelijks zichtbaar.

small %>%
  filter(kind == "SNV", pass) %>%
  mutate(set = ifelse(keep, "Gehouden", "Gefilterd")) %>% 
  dplyr::count(chrom, set) %>%
  mutate(
    per_mb = n / (chrom_len[chrom] / 1e6), 
    chrom = factor(chrom, levels = main_chromosomen)
  ) %>%
  ggplot(aes(x = chrom, y = per_mb, fill = set)) +
  geom_col(position = "stack", width = 0.7) +
  scale_fill_manual(values = c("Gehouden" = "#2b5c8f", "Gefilterd" = "#d95f02")) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major.x = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(color = "gray40", size = 10, margin = margin(b = 10)),
    legend.position = "bottom"
  ) +
  labs(
    x = NULL, 
    y = "SNVs per Mb", 
    fill = NULL, 
    title = "PASS SNVs per Mb per chromosoom",
    subtitle = "Een vlakker profiel is kenmerkend voor een goed tumor/controle-paar"
  )

Interpretatie: De verdeling is verre van vlak. chr2L heeft bijna 150 keer zoveel PASS-SNV’s per Mb als chr3R. Echte somatische mutaties horen niet op een paar chromosomen te hopen en andere over te slaan. Dit is meestal veroorzaakt door fouten bij batching, mismatched controle/tumor. In dit geval lijken de tumorsample en de normaal tot dezelfde batch te behoren (op basis van metadata), dus ik heb geen verklaring hiervoor.

\(\\\) \(\\\)

VAF grafiek: Op de x-as staat de VAF (variant allele frequency) in de tumor (van 0 tot 0,5), op de y-as het aantal PASS-SNVs. De rode stippellijn is de drempel van 0,2 zoals in de paper. Voor meer duidelijkheid over VAF zie (Smith et al. 2025).

snv <- small %>%
  filter(kind == "SNV", pass)

ggplot(snv, aes(t_vaf)) +
  geom_histogram(
    binwidth = 0.02,
    boundary = 0,
    fill = "grey40",
    colour = "white",
    linewidth = 0.2
  ) +
  geom_vline(
    xintercept = min_vaf,
    linetype = "dashed",
    colour = "red",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = min_vaf + 0.01,
    y = Inf,
    label = paste0("Minimum VAF = ", min_vaf),
    vjust = 1.5,
    hjust = 0,
    colour = "red"
  ) +
  scale_x_continuous(
    limits = c(0, .5),
    breaks = seq(0, 1, 0.1),
    expand = c(0, 0)
  ) +
  labs(
    x = "Variant Allel Fractie",
    y = "n PASS SNVs",
    title = "VAF per PASS SNV",
    subtitle = "Meeste varianten vallen beneden de 0.2 grens waarop gefilterd word."
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(colour = "grey40"),
    panel.grid.minor = element_blank()
  )

mean(snv$t_vaf < min_vaf, na.rm=T) * 100 # 64.9
## [1] 64.90773

Interpretatie: Een mutatie die in vrijwel alle cellen van een zuiver monster zit, heeft een VAF van ongeveer 0,5. Waarden van 0,1 tot 0,2 betekenen dat de variant maar in een klein deel van de cellen zit, of dat het om ruis of een mengsel van genotypen gaat. In deze data heeft 64.9% van de PASS SNVs een VAF van onder de 0.2. Ter vergelijking: in het artikel had ongeveer 93% van de varianten een allelfrequentie onder 0,2. Het percentage ligt hier lager, maar de richting is gelijk: de meeste calls hebben een lage VAF. De auteurs concluderen dat zulke varianten een tumor niet kunnen aansturen. Zonder schatting van de tumorpuriteit kunnen we niet onderscheiden of een lage VAF komt door een mutatie, een laag tumoraandeel, of door ruis en een verschil in achtergrond tussen tumor en controle.

\(\\\) \(\\\)

Structurele varianten

Wat de code doet: Voor Manta tellen we de ondersteunende readparen (PR) en split-reads (SR) in tumor en controle op. Voor Delly gebruiken we DV (readparen) + RV (junction reads). Een SV (structural variant) blijft wederom staan als hij PASS is, door de caller als somatisch gemarkeerd is, op een hoofdchromosoom ligt, minstens 5 ondersteunende reads in de tumor heeft (zoals in het artikel) en geen reads in de controle.

# Read structural variants from Manta or Delly
read_sv <- function(path, caller) {

  vcf <- readVcf(path, genome = "dm6")
  s <- who(vcf)
  g <- geno(vcf)
  inf <- info(vcf)

  if (caller == "Manta") {

    # Manta stores read-pair and split-read support as ref,alt
    get_support <- function(sample) {
      pr <- vapply(g$PR[, sample], \(x) if (length(x) >= 2) x[2] else NA, numeric(1))
      sr <- vapply(g$SR[, sample], \(x) if (length(x) >= 2) x[2] else NA, numeric(1))
      rowSums(cbind(pr, sr), na.rm = TRUE)
    }

    # Manta writes some variants twice, once for each side
    mate <- vapply(inf$MATEID, \(x) if (length(x)) x[1] else NA_character_, character(1))
    keep_pair <- is.na(mate) | names(vcf) < mate

    somatic <- TRUE

  } else {

    # Delly: variant read pairs + variant junction reads
    get_support <- function(sample) {
      g$DV[, sample] + g$RV[, sample]
    }

    keep_pair <- TRUE
    somatic <- inf$SOMATIC
  }

  tibble(
    caller = caller,
    type = as.character(inf$SVTYPE),
    chrom = as.character(seqnames(vcf)),
    pos = start(vcf),
    filter = as.character(rowRanges(vcf)$FILTER),
    somatic = somatic,
    t_support = get_support(s["tumor"]),
    n_support = get_support(s["normal"])
  )[keep_pair, ]
}

# Combine Manta and Delly calls
sv <- bind_rows(
  read_sv(manta_vcf, "Manta"),
  read_sv(delly_vcf, "Delly")
) %>%
  mutate(
    keep = filter == "PASS" &
           somatic &
           chrom %in% main_chromosomen &
           t_support >= sv_min_support &
           n_support == 0
  )

# summary van hoeveel blijven/ hoeveel droppen
sv %>%
  dplyr::count(caller, keep) %>%
  knitr::kable()
caller keep n
Delly FALSE 14
Delly TRUE 6
Manta FALSE 18
Manta TRUE 7 
sv %>%
  filter(keep) %>%
  arrange(chrom, pos) %>%
  knitr::kable()
caller type chrom pos filter somatic t_support n_support keep
Delly INS chr2L 323801 PASS TRUE 16 0 TRUE
Delly INV chr2L 1342705 PASS TRUE 7 0 TRUE
Manta DEL chr2L 11889472 PASS TRUE 8 0 TRUE
Delly DEL chr2R 3272989 PASS TRUE 11 0 TRUE
Delly DEL chr2R 5301902 PASS TRUE 12 0 TRUE
Manta DEL chr2R 5472475 PASS TRUE 9 0 TRUE
Manta DEL chr2R 5678814 PASS TRUE 7 0 TRUE
Manta BND chr2R 11790130 PASS TRUE 5 0 TRUE
Delly DEL chr2R 12303144 PASS TRUE 13 0 TRUE
Manta DEL chr2R 17387792 PASS TRUE 7 0 TRUE
Manta DEL chr3L 1318946 PASS TRUE 9 0 TRUE
Manta BND chr3L 24565325 PASS TRUE 7 0 TRUE
Delly DEL chrX 6263729 PASS TRUE 25 0 TRUE

Resultaat: De 13 gehouden SV’s: - Delly: insertie chr2L:323.801 (16 reads), inversie chr2L:1.342.705 (7 reads), deleties chr2R:3.272.989 (11), chr2R:5.301.902 (12), chr2R:12.303.144 (13) en chrX:6.263.729 (25).

  • Manta: deleties chr2L:11.889.472 (8), chr2R:5.472.475 (9), chr2R:5.678.814 (7), chr2R:17.387.792 (7), chr3L:1.318.946 (9) en twee translocaties (breekpunten bij chr2R:11.790.130 met 5 reads en chr3L:24.565.325 met 7 reads).

Interpretatie: In deze SVs zien we in de tumor ondersteunende reads en meerdere reads per SV. Dat is wat we van een somatische SV verwachten. Er zijn wel twee belangrijke punten van kritiek:

  1. Geen SNVs gevonden door beide callers.
  2. Sommige SVs worden ondersteund door maar weinig reads.

\(\\\) \(\\\)

Copy numbers

CNV grafiek: Elk punt is een segment. Op de x-as staat de positie binnen het chromosoom, op de y-as de log2-ratio uit CNVkit. Grijs is neutraal (rond de 0). De stippellijnen bij +0,585 en −0,6 zijn de drempels voor gain en loss uit het artikel (ratio 1,5 en ongeveer 0,66 berekend met logaritme). Een punt buiten de lijnen en met minstens 20 probes (zie config) wordt rood (gain) of blauw (loss).

cnv <- read.delim(cnvkit_cns) %>%
  filter(chromosome %in% main_chromosomen) %>%
  mutate(call = case_when(probes < cnv_min_probes ~ "neutral",
                          log2 >= cnv_gain ~ "gain",
                          log2 <= cnv_loss ~ "loss",
                          TRUE ~ "neutral"))

ggplot(cnv, aes(x = (start + end) / 2e6, y = log2, colour = call)) +
  geom_point(size = 1) +
  geom_hline(yintercept = c(cnv_loss, 0, cnv_gain), linetype = "dashed", colour = "orange4") +
  scale_colour_manual(values = c(gain = "red", loss = "blue", neutral = "grey55")) +
  facet_wrap(~ chromosome, scales = "free_x", ncol = 2) +
  labs(x = "Position (Mb)", y = "log2 ratio", title = "CNVkit segments", colour = NULL)

# cnv %>% dplyr::count(call) %>% knitr::kable()

Opmerking: De eerste versie van het cns call-bestand had een log2 die voor elk segment 0,212 hoger was dan in de segmentatie-uitvoer van CNVkit, waardoor neutrale regio’s op ongeveer +0,2 stonden en er vijf gains ontstonden. De stap cnvkit call is opnieuw uitgevoerd met –center-at -0.082 (een handmatige constante correctie die van alle log2-waarden wordt afgetrokken; in dit geval 0,082 erbij). Dat getal is de lengte-gewogen mediaan van de segmenten in de segmentatie-uitvoer. Daarna is de lengte-gewogen mediaan van het call-bestand precies 0, zoals bedoeld.

seg  <- read.delim("../../files/24-26/CNVkit_segs.cns")
call <- read.delim("../../files/24-26/CNVkit_CALL.cns") # oude call bestand

verschil <- call$log2 - seg$log2
summary(verschil)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2119  0.2119  0.2119  0.2119  0.2119  0.2119
seg <- read.delim("../../files/24-26/CNVkit_CALL.cns")

w <- seg$end - seg$start # lengte
o <- order(seg$log2) # log2 selecteren
baseline <- seg$log2[o][which(cumsum(w[o]) >= sum(w) / 2)[1]]   # gewogen mediaan
round(baseline, 3) # 0.13
## [1] 0.13
last <- baseline

( correctie <- round(mean(verschil - last),3) ) # correctie die handmatig moet ingevuld worden
## [1] 0.082

Resultaat: Van de 567 segmenten is er 1 als gain gecalld, 566 zijn neutraal en er is geen loss:

range(cnv$log2)
## [1] -0.538540  0.792039
cnv[which.min(cnv$log2), c("chromosome", "start", "end", "log2", "probes")]
##    chromosome    start      end     log2 probes
## 93       chrX 20290723 20310037 -0.53854     34
sum(cnv$probes < cnv_min_probes & (cnv$log2 >= cnv_gain | cnv$log2 <= cnv_loss))
## [1] 23

De log2-waarden lopen van −0,54 tot +0,79. Het diepste segment (chrX, 20,29 tot 20,31 Mb, log2 = −0,54) haalt de loss-drempel van −0,6 niet. Nog 23 segmenten liggen voorbij een drempel maar hebben minder dan 20 probes en tellen niet mee.

cnv %>%
  filter(call != "neutral") %>%
  mutate(start = round(start / 1e6, 2),
         end = round(end / 1e6, 2),
         log2 = round(log2, 3)) %>%
  dplyr::select(chromosome, start, end, log2, probes, call)
##   chromosome start   end  log2 probes call
## 1      chr3L 26.88 27.02 0.597     40 gain

De enigste ‘gain’ is in chr3L, 26,88 tot 27,02 Mb, log2 0,597, 40 probes.

Interpretatie: Er is geen duidelijke copynumbervariation. Het enige gecallde segment (chr3L) is een gain die maar net boven de drempel zit: 0,597 tegen een drempel van 0,585. Het is ook klein (ongeveer 0,14 Mb, 40 probes). Op zichzelf is dit segment niet onmogelijk, maar het zit zo dicht bij de grens dat een kleine verschuiving van de basislijn het in of uit de lijst haalt (wat wij in principe ook handmatig hebben gedaan). Het volgende segment (chr2L, 21,44 tot 21,54 Mb) komt op 0,556 en blijft onder de drempel. Dit komt overeen met het artikel, waar geen terugkerende kopiegetalveranderingen werden gevonden.

Samenvatting

De laatste tabel zet per tool het aantal PASS-calls naast het aantal uiteindelijk gehouden calls. “PASS” is wat de caller zelf heeft goedkeurd, “kept” is wat daadwerkelijk ook de extra filtering haalt. Voor CNVkit bestaat PASS niet (daarom NA) en is “kept” het aantal segmenten dat gecalld is.

bind_rows(
  small %>% 
    group_by(tool = paste("Strelka2", kind)) %>% 
    summarise(PASS = sum(pass), kept = sum(keep)),
  sv %>% 
    group_by(tool = paste(caller, "SV")) %>% 
    summarise(PASS = sum(filter == "PASS"), 
              kept = sum(keep)),
  tibble(tool = "CNVkit segments", PASS = NA, kept = sum(cnv$call != "neutral"))) %>%
  knitr::kable()
tool PASS kept
Strelka2 SNV 10025 65
Strelka2 indel 593 12
Delly SV 7 6
Manta SV 9 7
CNVkit segments NA 1

Conclusie en beperkingen

Conclusie: Er zijn candidate-somatische varianten in de tumor in vergelijking met de controle, maar geen daarvan is overtuigend bevestigd:

  • SVs: 13 kandidaten met 5 tot 25 ondersteunende reads in de tumor, maar geen overlapping met 2de caller.
  • SNVs en indels: na filteren blijven 65 SNV’s en 12 indels over. Met vooral calls op chromosomen chr2 en chrX lijken vooral verschillen tussen de 2 monsters aangetoond en niet mutaties.
  • Copynumber: na het rechtrekken van de basislijn blijft er 1 call over (een gain op chr3L, log2 = 0,597 bij een drempel van 0,585), en geen loss. Er zijn dus geen duidelijke kopiegetalveranderingen.

Qua resultaten sluit dit aan bij het artikel, dat geen terugkerende, hoogfrequente mutaties vond die de tumoren kunnen verklaren.

Om dan als laatst nog terug te komen op de onderzoeksvraag: Ja, we kunnen met deze pipeline constateren dat er wel degelijk variatie/afwijkingen bestaan in de behandelde toestand tegenover de controleconditie, maar deze kunnen niet direct gelinked worden aan de treatment (constant RNAi).

Beperkingen:

  1. Voor nu is alleen 1 sample-paar geanalyseerd. Met resultaten van het team is misschien meer te zeggen over het gehele experiment.
  2. De SRR accessions zijn niet goed geannoteerd door de originele auteurs; het koppelen van samples naar batch is gedaan op basis van metadata en order van de lijst, en lijken daarmee tot dezelfde batch te behoren.
  3. Tumorzuiverheid is niet berekend.
  4. Genannotatie is niet gedaan. (In dit geval niet nodig voor de conclusie, maar is interessant om de getroffen genen en hun funties te zien.)
  5. Gebruik van andere tools. Naast de vooraf besloten uitwisseling van MuTect naar Strelka2 zijn Breakdancer en CNVnator ook uitgewisseld voor modernere tools met zelfde functie. Deze veranderingen zijn zodanig groot in vergelijking met het artikel dat het gedane onderzoek niet meer een reproductie is.
  6. De eis voor geen alt-reads in de normaal voor een somatische variant is erg streng, de readdepth is in dit geval 109x en 1 alt read kan gemakkelijk noise zijn binnen de reads.

\(\\\) \(\\\)

Galaxy Workflow

Dit is een overview van de Galaxy Workflow bijbehorende de stappen van de pipeline en dit logboek.

(Als dit een wit blok is is dat jammer, de workflow is dan hier te vinden.)

\(\\\) \(\\\)

Jesse Postma, 507655, 09-10-2026, Genomics & Transcriptomics \(\\\)

Config & Versies

Tool versies

Tool Versie Waar gerund
fastp 0.24.0 & 1.3.7 BIN-netwerk (oud!) & Galaxy
BWA 0.7.18-r1243-dirty & 2.3+Galaxy BIN-netwerk & Galaxy
samtools 1.21 BIN-netwerk / NUC
Picard 3.1.1.1 BIN-netwerk (verzie onzeker) & Galaxy
Manta 1.6+galaxy9 Galaxy
Strelka2 2.9.10(+galaxy0) podman-container op NUC 412 & Galaxy
CNVkit 0.9.12+galaxy0 Galaxy
Delly 0.9.1+galaxy1 Galaxy
R sessionInfo 
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Debian GNU/Linux 13 (trixie)
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libmkl_rt.so;  LAPACK version 3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C           
##  [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C       
##  [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C        
## [10] LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 
## 
## time zone: Europe/Amsterdam
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggplot2_4.0.3               tidyr_1.3.2                
##  [3] dplyr_1.2.1                 VariantAnnotation_1.58.0   
##  [5] Rsamtools_2.28.0            Biostrings_2.80.2          
##  [7] XVector_0.52.0              SummarizedExperiment_1.42.0
##  [9] Biobase_2.72.0              GenomicRanges_1.64.0       
## [11] IRanges_2.46.0              S4Vectors_0.50.3           
## [13] Seqinfo_1.2.0               MatrixGenerics_1.24.0      
## [15] matrixStats_1.5.0           BiocGenerics_0.58.1        
## [17] generics_0.1.4              bslib_0.12.0               
## 
## loaded via a namespace (and not attached):
##  [1] KEGGREST_1.52.2          gtable_0.3.6             rjson_0.2.23            
##  [4] xfun_0.61                lattice_0.23-1           vctrs_0.7.3             
##  [7] tools_4.6.1              bitops_1.1-0             curl_8.0.0              
## [10] parallel_4.6.1           tibble_3.3.1             AnnotationDbi_1.74.0    
## [13] RSQLite_3.53.3           blob_1.3.0               pkgconfig_2.0.3         
## [16] Matrix_1.7-6             BSgenome_1.80.0          RColorBrewer_1.1-3      
## [19] S7_0.2.2                 cigarillo_1.2.1          lifecycle_1.0.5         
## [22] farver_2.1.2             compiler_4.6.1           codetools_0.2-20        
## [25] htmltools_0.5.9          sass_0.4.10              RCurl_1.98-1.20         
## [28] yaml_2.3.12              pillar_1.11.1            crayon_1.5.3            
## [31] jquerylib_0.1.4          BiocParallel_1.46.0      DelayedArray_0.38.2     
## [34] cachem_1.1.0             abind_1.4-8              tidyselect_1.2.1        
## [37] digest_0.6.39            purrr_1.2.2              restfulr_0.0.17         
## [40] labeling_0.4.3           fastmap_1.2.0            grid_4.6.1              
## [43] cli_3.6.6                SparseArray_1.12.2       magrittr_2.0.5          
## [46] S4Arrays_1.12.0          GenomicFeatures_1.64.0   utf8_1.2.6              
## [49] dichromat_2.0-1          XML_3.99-0.24            withr_3.0.3             
## [52] scales_1.4.0             bit64_4.8.6              rmarkdown_2.32          
## [55] httr_1.4.9               bit_4.6.0                otel_0.2.0              
## [58] png_0.1-9                memoise_2.0.1            evaluate_1.0.5          
## [61] knitr_1.52               BiocIO_1.22.0            rtracklayer_1.72.0      
## [64] rlang_1.3.0              glue_1.8.1               DBI_1.3.0               
## [67] rstudioapi_0.19.0        jsonlite_2.0.0           R6_2.6.1                
## [70] GenomicAlignments_1.48.0 fs_2.1.0

Referenties

Cai, Lei, Wei Yuan, Zhou Zhang, Lin He, and Kuo Chen Chou. 2016. “In-Depth Comparison of Somatic Point Mutation Callers Based on Different Tumor Next-Generation Sequencing Depth Data.” Scientific Reports 2016 6:1 6 (November): 36540–40. https://doi.org/10.1038/srep36540.
GeneAnalyses Editorial Team. 2026. “Variant Allele Fraction: What the Percentage on a Report Means.” GeneAnalyses, July. https://www.geneanalyses.com/articles/variant-allele-fraction-explained.
Kim, Sangtae, Konrad Scheffler, Aaron L. Halpern, et al. 2018. “Strelka2: fast and accurate calling of germline and somatic variants.” Nat. Methods 15 (August): 591–94. https://doi.org/10.1038/s41592-018-0051-x.
Koboldt, Daniel C., Qunyuan Zhang, David E. Larson, et al. 2012. “VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing.” Genome Res., ahead of print, February. https://doi.org/10.1101/gr.129684.111.
Parreno, V., V. Loubiere, B. Schuettengruber, et al. 2024. “Transient loss of Polycomb components induces an epigenetic cancer fate.” Nature 629 (May): 688–96. https://doi.org/10.1038/s41586-024-07328-w.
Smith, Adam C., Hubert Tsui, Sila Usta, and Jose-Mario Capo-Chichi. 2025. “What the VAF? A guide to the interpretation of variant allele fraction, percent mosaicism, and copy number in cancer.” Mol. Cytogenet. 18 (1): 13. https://doi.org/10.1186/s13039-025-00718-3.