Transcriptomische analyse gebaseerd op bulk RNA-seq gegevens

Summary

Het huidige protocol stelt een volledige pijplijn vast voor het analyseren van het proces van bulk RNA-seq van ruwe data tot functionele verrijkingsanalyse.

Abstract

Niet-alcoholische vettige lever (NAFL) wordt meestal beschouwd als een goedaardige aandoening; echter, zodra het overgaat naar niet-alcoholische steatohepatitis (NASH), lopen patiënten een aanzienlijk verhoogd risico op het ontwikkelen van leverziekte in het eindstadium. Veel studies proberen het moleculaire mechanisme achter de overgang van NAFL naar NASH te verduidelijken. High-throughput sequencingtechnologieën (zoals bulk RNA-seq) hebben onderzoekers een dieper inzicht gegeven door het transcriptoom te onderzoeken, wat de expressie van moleculen, activatie van signaalroutes en andere factoren die samenhangen met ziekteprogressie onthult. Er is een schat aan open source data beschikbaar voor onderzoekers om te analyseren om potentiële doelwitten voor ziektebehandeling te identificeren. Gerelateerd onderzoek wordt echter beperkt door het ontbreken van een efficiënt en betrouwbaar proces voor upstream-analyse van het transcriptoom. Hier wordt een zeer reproduceerbare en gebruiksvriendelijke upstream-analyse en daaropvolgende gerelateerde differentiële genanalysepijplijn aangeboden om gestandaardiseerde verwerking en diepgaande parsing van private of publieke data te bereiken. De pijplijn is verdeeld in vier stappen: (1) kwaliteitscontrole van data; (2) genmapping; (3) differentiële genanalyse; en (4) functionaalanalyse. Dit proces heeft als doel de moleculaire mechanismen van ziektetransformatie te onthullen en onderzoekers te helpen potentiële geneesmiddeldoelen en therapeutische benaderingen te screenen via de analyse van Bulk RNA-seq data.

Introduction

Niet-alcoholische leververvetting (NAFLD) is wereldwijd de meest voorkomende chronische leverziekte en treft meer dan een kwart van de bevolking. De incidentie is de afgelopen decennia dramatisch toegenomen met ^1,2,3. De groeiende ziektelast, vooral de meer gevorderde vorm, niet-alcoholische steatohepatitis (NASH), vormt een grote wereldwijde gezondheidsuitdaging en een zware economische last⁴. De eerste fase van NAFLD is niet-alcoholische levervetting (NAFL), die gepaard gaat met ontstekingen en fibrose die kunnen overgaan in NASH. Deze laatste verhoogt aanzienlijk het risico op progressie naar leverziekte in het eindstadium, waaronder cirrose en hepatocellulair carcinoom (HCC)^5,6,7. HCC-incidentie en mortaliteit zijn geassocieerd met een toename van NASH ^8,9, en het wordt verwacht dat NAFLD/NASH tegen 2030 de belangrijkste indicatie voor levertransplantatie zal^worden. De klinische progressie van NAFLD is echter zeer heterogeen¹¹, wat de ontwikkeling van relevante geneesmiddelen¹² ernstig belemmert, waardoor het bijzonder belangrijk is om de betrokken moleculaire mechanismen nauwkeurig te onderzoeken.

Bulk-RNA-seq-gebaseerde verwerving van cellulaire samenstellingsinformatie kan de pathogenese van verschillende ziekten aanzienlijk verduidelijken. In de afgelopen decennia zijn talrijke bulk-RNA-seq-studies uitgevoerd bij modelorganismen en mensen om verschillen in genexpressie in NASH-progressie 13,14,15 te verduidelijken en zo nieuwe therapeutische doelwitten voor interventie te identificeren. Op basis van bulk RNA-seq-analyse ontdekten Xiong et al. dat niet-parenchymale cellen (NPC's) in de lever betrokken zijn bij processen zoals de vorming van extracellulaire matrixs en celadhesie, wat bijdraagt aan de progressie van NASH¹⁶. Li et al. toonden aan dat het hepatische Wilms' tumor 1-associatieve eiwit (WTAP) in hepatocyten de ophoping en ontsteking van ectopische lipiden reguleert, waardoor NASH-vorming¹⁷ wordt bevorderd. Hoewel bulk RNA-seq-analyse een krachtig hulpmiddel is om de mechanismen van NASH te verduidelijken, zijn de resultaten zeer gevoelig voor de kwaliteit van upstream-gegevens. De heterogeniteit van experimentele operaties en analyseprocessen stroomopwaarts kan de betrouwbaarheid van data ernstig aantasten, waardoor echte biologische informatie wordt verhuld en de nauwkeurigheid van latere analyses wordt belemmerd. Daarom is het belangrijk om een set gestandaardiseerde upstream analyseprocedures op te stellen.

In vergelijking met single-cell RNA-sequencing (scRNA-seq) biedt bulk RNA-seq verschillende duidelijke voordelen in zowel experimenteel ontwerp als praktische toepassingen. Hoewel scRNA-seq de identificatie van cellulaire heterogeniteit op enkelcelniveau mogelijk maakt en nauwkeurige analyse van celtype-specifieke transcriptiekenmerken mogelijk maakt, gaat het gepaard met hoge kosten, complexe gegevensverwerkingsvereisten en beperkte gevoeligheid voor het detecteren van laag-abundantie transcripten¹⁸. Daarentegen biedt bulk RNA-seq een hogere sequencingdiepte, lagere kosten en een grotere steekproefdoorvoer, waardoor het bijzonder geschikt is voor populatieniveau-differentiële genexpressie-analyses en de verkenning van moleculaire mechanismen¹⁹. Daarom blijft bulk RNA-seq, wanneer het wordt geleid door gestandaardiseerde analytische workflows, een efficiënte, kosteneffectieve en robuuste benadering voor het onderzoeken van de moleculaire basis van complexe ziekten.

Dit protocol is specifiek ontworpen voor bulk RNA-seq datasets afkomstig van menselijke weefsels met een hoge RNA-integriteit (RIN ≥ 7,0) en voldoende input-RNA (≥ 500 ng per monster). Om een betrouwbare uitvoering van uitlijnings- en kwantificatiestappen te waarborgen, wordt een lokaal werkstation aanbevolen met ten minste een 10-core CPU, 32 GB RAM en minimaal 200 GB vrije schijfruimte. Voortbouwend op deze eisen biedt het protocol een efficiënte en gebruiksvriendelijke analytische workflow, inclusief gedetailleerde operationele instructies en gestandaardiseerde parameterconfiguraties, om te voldoen aan de behoeften van onderzoekers die grootschalige transcriptomische data analyseren.

Protocol

Voor demonstratiedoeleinden werd de publiek beschikbare dataset PRJNA1023502 gegenereerd door Lan Bai et al. gebruikt om elke stap van zowel upstream- als downstream-analyses te illustreren²⁰. Omdat deze dataset afkomstig is uit de open-access NCBI SRA-database, zijn geen extra toestemmingen of ethische goedkeuringen vereist. Zie de Materiaaltabel om alle vereiste software- en R-pakketversies te verifiëren. De publiek beschikbare dataset bestaat PRJNA1023502 uit 6 niet-NASH, 6 NAFL- en 6 NASH-lever-RNA-seq monsters. In dit protocol werd de dataset gebruikt om alle stappen van de bulk RNA-seq workflow te demonstreren, inclusief gegevensopvraging uit de SRA-database, kwaliteitscontrole (fastp), uitlijning (HISAT2), kwantificatie (featureCounts) en downstream differentiële expressie- en functionele verrijkingsanalyses.

1. SRA-toolkitinstallatie

1. Bezoek de officiële website van de SRA Toolkit en download versie 3.2.1.

2. Openbare datadownload

1. Verkrijg het SRA-nummer.
   1. Toegang hiernaar in artikelsupplement²⁰, de sectie Data Availability, of door te zoeken naar trefwoorden in de NCBI SRA-database.
   2. Typ prefetch in de terminal om het te downloaden.

3. Generatie van gentellingsmatrix

1. Definieer paden en instellingen zoals hieronder beschreven:
   REFERENCE=~/reference/human/GRCh38/GRCh38.primary_assembly.genome.fa
GTF=~/reference/human/GRCh38/gencode.v44.annotation.gtf
INDEX=~/reference/human/GRCh38/GRCh38_index
FASTQ_DIR=~/SRA_tutorial/fastq
OUT_FASTP=~/RNAseq/fastp
OUT_HISAT2=~/RNAseq/hisat2
OUT_COUNTS=~/RNAseq/counts
​mkdir -p $FASTQ_DIR $OUT_FASTP $OUT_HISAT2 $OUT_COUNTS
   1. Download het menselijke referentiegenoom (GRCh38 primaire assemblage) en het bijbehorende genannotatiebestand (gencode.v44, uitgebreide genannotatieset) uit de officiële GENCODE-database (https://www.gencodegenes.org/human/).
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
​for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
2. Hernoem SRA-bestanden.
   1. Zorg ervoor dat alle SRA-bestanden de .sra-bestandsextensie gebruiken om nauwkeurige herkenning en verwerking door downstream-tools te vergemakkelijken.
   2. Voer de volgende commando's uit:
      for f in SRR*; do [[ ! $f =~ \.sra$ ]] && mv "$f" "$f.sra"; done
3. Voer het volgende commando uit om SRA-bestanden om te zetten naar FASTQ-formaat:
   for f in *.sra; do fasterq-dump "$f" --split-files -O $FASTQ_DIR - e 20; done
4. Bouw een HISAT2-index op basis van het referentiegenoom.
   hisat2-build $REFERENCE $INDEX
5. Verwerk FASTQ-bestanden met visuele voortgangsmonitoring.
   for fq in $FASTQ_DIR/*.fastq; do
​sample=$(basename "$fq" .fastq)
   1. Voor paired-end sequencinggegevens voer je het volgende commando uit
      for fq1 in $FASTQ_DIR/*_1.fastq; do
sample=$(basename "$fq1" _1.fastq)
​fq2=$FASTQ_DIR/${sample}_2.fastq
   2. Voer kwaliteitscontrole en filtering uit met FASTP en voer de resultaten uit naar de $OUT_FASTP-map. Voor single-end sequencingdata voer je het volgende commando uit:
      fastp \
-i "${fq}" \
-o $OUT_FASTP/${sample}.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
-w 20
      Voor paired-end sequencinggegevens voer je het volgende commando uit:
      fastp \
-i "${fq}" \ -I "$fq2" \
-o $OUT_FASTP/${sample}_1.clean.fastq \
-O $OUT_FASTP/${sample}_2.clean.fastq \
-h $OUT_FASTP/${sample}.html \
-j $OUT_FASTP/${sample}.json \
​-w 20
   3. Voor elk monster inspecteert u het HTML-rapport op visuele kwaliteitscontrole en raadpleegt u het JSON-rapport voor gestructureerde metrics, waaronder leeskwaliteit, GC-distributie, leesredundantie en adapterverontreiniging. Zorg ervoor dat alle monsters voldoen aan de volgende kwaliteitsdrempels: Q30 ≥ 85%, adapterverontreinigingsgraad < 5% en geen afwijkende GC-patronen. Sla alle uitvoerbestanden op in de $OUT_FASTP-map.
   4. Voer sequence alignment uit met HISAT2, genereer SAM-bestanden en voer deze uit naar de $OUT_HISAT2-map. Voor single-end sequencingdata voer je het volgende commando uit:
      hisat2 -p 20 \ -x $INDEX \-U $OUT_FASTP/${sample}.clean.fastq \
-S $OUT_HISAT2/${sample}.sam
      Voor paired-end sequencinggegevens voer je het volgende commando uit:
      hisat2 -p 20 \-x $INDEX \-1 $OUT_FASTP/${sample}_1.clean.fastq \
-2 $OUT_FASTP/${sample}_2.clean.fastq \
​-S $OUT_HISAT2/${sample}.sam
   5. Voer HISAT2-uitlijning uit van de kwaliteitsgefilterde leesopdrachten naar het menselijke referentiegenoom, met behulp van een index die is opgebouwd uit het GRCh38 primaire assemblagebestand FASTA. Geef voorbeeldcommando's voor zowel single-end als paired-end sequencingdata.
   6. Zet SAM om naar BAM, sorteer en indexeer.
      samtools view -@ 20 -bS $OUT_HISAT2/${sample}.sam \
| samtools sort -@ 20 -o $OUT_HISAT2/${sample}.sorted.bam
samtools index $OUT_HISAT2/${sample}.sorted.bam
​done
   7. Genereer een SAM-bestand voor elk voorbeeld en zet dit om naar een gesorteerd en geïndexeerd BAM-bestand met SAMtools. Voor hoogwaardige menselijke transcriptomische monsters (RNA-integriteitsgetal, RIN ≥ 8,0) moet de totale uitlijningssnelheid meer dan 85% worden gezorgd; voor typische bulk RNA-seq-monsters (RIN ≥ 7,0) wordt de uitlijningspercentages ≥ 70% acceptabel beschouwd.
   8. Voer genkwantificatie uit met featureCount.
      featureCounts -T 20 -p -s 0 \
-a $GTF \
-o $OUT_COUNTS /${sample}.counts.txt \
$OUT_HISAT2/${sample}.sorted.bam
Done
   9. Controleer het tab-gescheiden uitvoerbestand (*.counts.txt) en het samenvattende rapport (*.counts.txt.summary) dat door featureCounts wordt gegenereerd voor elk voorbeeld. Zorg ervoor dat de leestoewijzingssnelheid de gebruikelijke drempel van ≥70% voor human bulk RNA-seq haalt; Een aanzienlijk lagere snelheid kan wijzen op een niet-matchende strengoriëntatie, annotatieproblemen of slechte uitlijningskwaliteit. Gebruik de parameter -s 0 voor deze niet-strengspecifieke RNA-seq dataset. Voor streng-specifieke bibliotheken vervang -s 0 door -s 1 of -s 2 in het commando.
6. Genereer de gentellingsmatrix zoals hieronder beschreven.
   1. Initialiseer de countsmatrix met gen-ID's en tellingen van het eerste monster.
      cut -f1 $(ls $OUT_COUNTS/*.counts.txt | head -1) > all_counts.txt
   2. Voeg iteratief tellingen van elke sample op in één matrix.
      for f in $OUT_COUNTS/*.counts.txt; do
cut -f7 "$f" | paste all_counts.txt - > tmp && mv tmp
all_counts.txt
​done
   3. Voeg een headerregel met voorbeeld-ID's toe aan de countsmatrix.
      samples=$(ls *.counts.txt | sed 's/.counts.txt//' | paste -sd "\t")
echo -e "Geneid\t$samples" | cat - all_counts.txt > counts_matrix.txt
   4. Haal genlengtes uit het GTF-bestand (som van exonlengtes per gen).
      awk '$3=="exon"{match($0,/gene_id "([^"]+)"/,a); if(a[1]!=""){len=$5-$4+1; gene_len[a[1]]+=len}} END{print "GENE_ID\tLENGTH"; for(g in gene_len) print g"\t"gene_len[g]}' \$GTF > gene_length.txt

4. Verwerking van ruwe telmatrixen en genannotatie

1. R-omgevingsinitialisatie en het laden van het relevante R-pakket.
   1. Gebruik de install.packages()-functie om tidyverse, ggplot2, ggrepel, RColorBrewer, ggridges, FactoMineR te installeren. Het tidyverse-pakket wordt gebruikt voor databewerking en plotting; Het ggplot2-pakket wordt gebruikt voor visualisatie; Het Ggrepel-pakket wordt gebruikt voor niet-overlappende tekstlabels; het RColorBrewer-pakket wordt gebruikt voor kleurpaletten; Het GGRIDGES-pakket wordt gebruikt voor Ridge-percelen; het FactoMineR-pakket wordt gebruikt voor PCA- en multivariate analyse.
   2. Gebruik de functie BiocManager::install() om biomaRt, DESeq2, clusterProfiler te installeren. Het biomaRt-pakket wordt gebruikt voor genannotatie vanuit Ensembl; het DESeq2-pakket wordt gebruikt voor differentiële expressieanalyse; het clusterProfiler-pakket wordt gebruikt voor functionele verrijkingsanalyse.
   3. Lees de originele telmatrix die in stap 2 is gegenereerd, gebruik biomaRt om Ensembl-ID te koppelen aan de HGNC-gennaam voor latere analyse, verwijder de dubbele Geneid-kolom en genereer de gereinigde telmatrix. Gebruik de originele countsmatrix (counts_matrix.csv) als input, met Ensembl-gen-ID's als rijen en monsters als kolommen. Bewaar de gereinigde tellingmatrix als output (clean_counts_SRA.csv), met HGNC-gensymbolen als rijen en monsters als kolommen.
      mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
id_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = exprSet$GeneID,
mart = mart)
exprSet <- exprSet %>%
left_join(id_map, by = c("GeneID" = "ensembl_gene_id")) %>%
filter(!is.na(hgnc_symbol), hgnc_symbol != "") %>%
distinct(hgnc_symbol, .keep_all = TRUE) %>%
column_to_rownames("hgnc_symbol")

5. Kwantificatie van genexpressie

OPMERKING: Raadpleeg het aanvullende dossier 1 voor het gedetailleerde script.

1. Voer het volgende commando uit om fragmenten per kilobase transcript per miljoen mapped reads (FPKM) en transcripties per miljoen (TPM) te berekenen.
   counts <- read.csv("output/clean_counts_SRA.csv", header=TRUE, row.names=1)
gene_len <- read.delim("data/gene_length.txt", header=FALSE, col.names=c("gene_symbol","length"))
gene_len <- gene_len %>% distinct(gene_symbol, .keep_all=TRUE)
rownames(gene_len) <- gene_len$gene_symbol
gene_len <- gene_len[match(rownames(counts), gene_len$gene_symbol),]
length_bp <- gene_len$length
fpkm <- (counts / length_bp) * 1e9 / colSums(counts)
write.csv(fpkm, "output/clean_fpkm_SRA.csv")
tpm <- (counts / length_bp) / colSums(counts / length_bp) * 1e6
write.csv(tpm, "output/clean_tpm_SRA.csv")

6. Voorbeeldclustering en verschilvisualisatie

1. Beoordeel de datakwaliteit door monsterclustering te onderzoeken via PCA. Raadpleeg Supplementair Dossier 1 voor het gedetailleerde script.
2. Voor PCA-visualisatie gebruik de onderstaande code.
   gene.pca <- PCA(exprSet, ncp = 2, scale.unit = TRUE, graph = FALSE)
ggplot(pca_sample, aes(x = Dim.1, y = Dim.2)) +
geom_point(aes(color = group)) +
labs(x = paste('PC1:', pca_eig1, '%'),
y = paste('PC2:', pca_eig2, '%'))

7. Analyse van differentiële expressie en visualisatie van resultaten

OPMERKING: Raadpleeg het aanvullende dossier 1 voor het gedetailleerde script.

1. Voer de volgende commando's uit om de DESeq2-dataset te construeren en stel de drempel in. Alleen genen met totale tellingen > 1 over alle monsters werden behouden vóór de modellering. Identificeer differentieel expressieve genen (DEGs) met behulp van de significantiedrempels van aangepaste p-waarde < 0,05 en |log2FC| > 0,5.
   dds <- DESeq(DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~group)); sizeFactors(dds); res <- results(dds); dds <- dds[rowSums(counts(dds)) > 1,]
dd1 <- results(dds, contrast = contrast, alpha = 0.05)
dd2 <- lfcShrink(dds, contrast = contrast, res = dd1, type = "ashr")
2. Visualiseer differentieel expressieve genen met een vulkaangrafiek.
   ggplot(data = data, aes(x = log2FoldChange, y = -log10(padj))) +
geom_point(aes(color = group), alpha = 1, size = 1.2) +
geom_hline(yintercept = -log10(0.05), lty = 4) +
geom_vline(xintercept = c(-0.5, 0.5), lty = 4) +
geom_text_repel(data = subset(data, abs(log2FoldChange) >= 1.5 & padj < 0.05),
aes(label = gene_id))

8. Voer functionele verrijkingsanalyse en visualisatie uit

OPMERKING: Raadpleeg het aanvullende dossier 1 voor het gedetailleerde script.

1. Visualiseer de KEGG-verrijkingsresultaten.
   EGG <- enrichKEGG(gene = gene$ENTREZID, organism = 'hsa',
pvalueCutoff = 0.05, qvalueCutoff = 0.05)
ggplot(symboldata, aes(richFactor, Description)) +
geom_point(aes(color = p.adjust, size = Count))
2. Visualiseer de resultaten van GO-verrijking. Pas de visualisatie aan door de boogresolutie (n=500), kleurafbeeldingen en facetconfiguratie aan te passen volgens experimentele specificaties.
   ego <- enrichGO(gene = gene$ENTREZID, OrgDb = "org.Hs.eg.db", ont = "ALL",
pvalueCutoff = 0.05, qvalueCutoff = 0.05, pAdjustMethod = "BH")
ggplot(df) +
ggforce::geom_link(aes(x = 0, y = Description, xend = -log10(p.adjust),
yend = Description, color = ONTOLOGY), n = 500, show.legend = FALSE) +
facet_wrap(~ONTOLOGY, scales = "free", ncol = 1)
3. Visualiseer de resultaten van gensetverrijkingsanalyse (GSEA). Pas het aantal paden, dichtheidsrugschaal en laagtransparantie aan om verschillende effectgrootteverdelingen te accommoderen.
   genelist <- sort(res$log2FoldChange, decreasing = TRUE)
names(genelist) <- rownames(res)
hallmarks <- read.gmt('resource/h.all.v2023.2.Hs.symbols.gmt')
y <- GSEA(genelist, TERM2GENE = hallmarks, pvalueCutoff = 0.05)
gsearesult <- yd %>% arrange(desc(NES)) %>% slice_head(n = 10)
ggplot(gsearesult, aes(x = logFC, y = Description, fill = -log10(pvalue))) +
geom_density_ridges(alpha = 0.8, scale = 0.8) +
geom_point(aes(size = abs(NES), x = -0.4, color = NES)) +
scale_fill_distiller(palette = 'Spectral') +
scale_color_distiller(palette = 'Reds') +
scale_size_continuous(range = c(2, 6))

Results

De upstream analyseworkflow voor bulk RNA-seq wordt geïllustreerd in Figuur 1A. Deze workflow voert sequentieel de volgende belangrijke stappen uit op een Linux-platform: ten eerste wordt rigoureuze kwaliteitscontrole van ruwe sequencinggegevens uitgevoerd met behulp van fastp om laagwaardige reads en adaptersequenties te verwijderen; vervolgens lijnt HISAT2 hoogwaardige leesopdrachten uit op het referentiegenoom, waarbij Samtools de uitlijningsbestanden converteert en sorteert; ten slotte v...

Discussion

Bulk RNA-seq data-analyse wordt gekenmerkt als een interdisciplinaire taak die genomica, bio-informatica, statistiek en informatica integreert. Een volledige analytische workflow omvat meerdere upstream- en downstreamstappen, waaronder preprocessing van ruwe data, kwaliteitscontrole, sequentie-uitlijning, kwantificatie op genniveau, datanormalisatie, differentiële expressie-analyse en biologische interpretatie. Van deze stappen is het nauwkeurig omzetten van ruwe sequencing-reads naar ee...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
biomaRt	Bioconductor	2.64.0	Genoomannotatie van Ensembl
clusterProfiler	Bioconductor	4.16.0	Functionele verrijkingsanalyse
DESeq2	Bioconductor	1.48.1	Differentiële expressieanalyse
FactoMineR	AgroParisTech	2.11.0	PCA en multivariate analyse
fastp	OpenGene	1.0.1	Controle en filtering van FASTQ-gegevens
FeatureCounts	Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research	2.0.0	Tell het aantal reads dat op elk gen is gemapped voor genexpressie kwantificering
ggplot2	Posit	3.5.2	Datavisualisatie
ggrepel	Kamil Slowikowski	0.9.6	Niet-overlappende tekstlabels
ggridges	Claus O. Wilke	0.5.6	Maak ridgeline-plots
HISAT2	Johns Hopkins University	2.2.1	Align de gefilterde hoge kwaliteit reads op het referentiegenoom
R	R Core Team	4.5.0	Een omgeving voor gegevensberekening, analyse en visualisatie
RColorBrewer	Erich Neuwirth	1.1.3	Kleurenpaletten voor plotten
samtools	Large Scale Genomics work stream	1.22.0	Converteren en verwerken van SAM-bestanden voor efficiënte opvraging en toegang
SRA Toolkit	National Center for Biotechnology Information	3.2.1	Verkrijgen en voorverwerken van ruwe sequencinggegevens van de NCBI SRA-database