Skapa ett eget nukleotidfrekvensdiagram
Nu är det dags att titta närmare på nukleotidfrekvensen per cykel. Det bästa sättet är att visualisera datan. Vanligtvis är de första cyklerna lite slumpmässiga, varefter nukleotidfrekvenserna stabiliserar sig.
Den här övningen använder den fullständiga fastq-filen SRR1971253 med viss förbearbetning gjord åt dig:
library(ShortRead)
fqsample <- readFastq(dirPath = "data",
pattern = "SRR1971253.fastq")
# extract reads
abc <- alphabetByCycle(sread(fqsample))
# Transpose nucleotides A, C, G, T per column
nucByCycle <- t(abc[1:4,])
# Tidy dataset
nucByCycle <- nucByCycle %>%
as_tibble() %>% # convert to tibble
mutate(cycle = 1:50) # add cycle numbers
Din uppgift är att skapa ett diagram över nukleotidfrekvens per cykel med hjälp av funktioner från tidyverse!
Den här övningen är en del av kursen
Introduktion till Bioconductor i R
Övningsinstruktioner
- Använd
glimpse()på objektetnucByCycleför att få en överblick över datan. - Pivotera nukleotidbokstäverna i
alphabetmedpivot_longer()och skapa en ny kolumncount. - Skapa ett linjediagram med
cyclepå x-axeln ochcountpå y-axeln, färgat efteralphabet.
Interaktiv övning med praktiskt arbete
Testa den här övningen genom att slutföra den här exempelkoden.
# Glimpse nucByCycle
___
# Create a line plot of cycle vs. count
nucByCycle %>%
# Gather the nucleotide letters in alphabet and get a new count column
pivot_longer(-cycle, names_to = ___, values_to = ___) %>%
ggplot(aes(x = ___, y = ___, color = ___)) +
geom_line(size = 0.5 ) +
labs(y = "Frequency") +
theme_bw() +
theme(panel.grid.major.x = element_blank())