1 Databaser

1.1 Litt bakgrunn

Ta først en titt på ukens forelesning. Det er spesielt hvordan vi søker opp data hos NCBI som er viktige for oppgavene nedenfor:

Hos NCBI finner du også massevis av dokumenter og filmer som viser hvordan du kan gå fram. Se for eksempel siden https://www.ncbi.nlm.nih.gov/home/learn/. Eksempel på nyttige filmer, med linker til annet på YouTube:

1.2 Pondus

Gjør et globalt søk på bilbane-genet hos NCBI. Hvor mange treff får du?

1.3 Nitrogen-omsetning i jord

I forbindelse med studier av nitrogen-omsetning i jord vil vi bygge opp en oversikt over hvilke varianter som finnes av et spesielt gen, og i hvilke bakterier de forekommer. Vi vil her fokusere på et gen som på engelsk kalles nitrite reductase og som i tillegg har et kort-navn nirK. Dette er et kodende gen. Søk med denne informasjonen i databasen Protein hos NCBI. Hvor mange treff får du?

Trolig er mange treff fra bakterier som er delvis ukjente, angitt som uncultured et-eller-annet. Prøv å utvide søke-teksten slik at du unngår treff mot alt som er uncultured. HINT: Bruk den logiske operatoren NOT.

Vi er ikke interessert i treff hos veldig korte deler av proteinet. Legg inn som et krav i søke-teksten at lengden skal være mellom 200 og 2000. HINT: teksten 200:2000[slen] angir at vi vil ha sekvenslengden i intervallet 200 til 2000.

Last ned alle treff som en FASTA-fil og lagre den på din PC. Vi skal bruke denne filen senere.



2 Filmer om sekvensering

Det finnes mange filmsnutter på nettet som illustrerer hvordan de ulike sekvenserings-teknologiene fungerer. Vi har valgt å fokusere på Illumina-teknologien på forelesningene, og her er en link du kan starte med:

Illumina: https://www.youtube.com/watch?v=fCd6B5HRaZ8

Du vil lett finne flere filmer på YouTube med liknende innhold, bruk litt tid på å studere disse. Illumina-teknologien er den mest utbredte akkurat nå. Den er ganske billig, og gir veldig presise sekvenser, det vi kaller reads. Derimot er disse read’ene ganske korte, vanligvis ikke mer enn 300 baser lange.

Det er selvsagt andre teknologier også, og typisk for de nyeste er at de lager mye lengre reads (tusener baser). Dette er en stor fordel i mange sammenhenger. Dessverre er de mer upresise enn Illumina, men dette endrer seg allerede mye. Her er de to vanlige teknologiene:

Nanopore: https://www.youtube.com/watch?v=E9-Rm5AoZGw

PacBio: https://www.youtube.com/watch?v=v8p4ph2MAvI



3 Inspisere data i FastQC

Data direkte fra sekvenserings-maskinen kommer vanligvis på det som kalles fastq formatet. Slike filer har da typisk endelsen .fastq. Fra Illumina får vi alltid et par av filer ettersom alle fragmenter sekvenseres fra begge ender, det vi kaller paired-end sekvensering. Selv om .fastq filer er tekstfiler, så er de ofte store og litt for “rotete” til å inspisere manuelt. Til det kan vi heller bruke en programvare som FastQC.

Ta først en kikk på videoen som viser hvordan du installerer programvaren FastQC:

Dersom du må installere Java for å få dette til å fungere, så følg anvisningene i fila INSTALL.txt som er i arkivet du pakker opp.

3.1 Hel-genom sekvensering av bakterie

La oss inspisere noen .fastq filer. Last ned to filer med data fra en Illumina MiSeq shotgun sekvensering av et bakterie-genom (høyre-klikk, Lagre…):

GenomeShotgun_R1.fastq.gz

GenomeShotgun_R2.fastq.gz

Det er lurt å lage en mappe med navn BIN210, et sted nær toppen av fil-treet, og legge filene i denne, eventuelt i en under-mappe i denne (f.eks. med navn data). Dette er paired-end data, og de opptrer alltid som et fil-par. Fila med R1 i navnet er reads fra runde 1, og R2 fra runde 2, se filmene over. Åpne filene fra FastQC, her er en filmsnutt som viser hvordan du går fram:

Husk at dette er ‘shotgun-sekvensering’ av et bakterie-genom. Det betyr at de ulike sekvensene er biter fra tilfeldige steder rundt omkring på genomet. Disse dataene er av gjennomgående ganske høy kvalitet, men du vil likevel finne noen forskjeller, for eksempel mellom R1-reads og R2-reads. Hva er ulikt? Hvorfor?

3.2 Amplicon-sekvensering av et og samme gen fra mange bakterier

Last ned filene

Amplicon16S_R1.fastq.gz

Amplicon16S_R2.fastq.gz

Åpne de i FastQC på samme måte som over.

I dette forsøket har vi sekvensert ett eneste gen, nemlig 16S-genet, hos en mengde ukjente bakterier. Vi har kopiert opp en spesiell bit av dette genet fra alle bakterier i en prøve fra et miljø, ved bruk av PCR og med primere som matcher akkurat dette genet. Deretter er de amplifiserte fragmentene sekvensert med Illumina MiSeq. Det betyr at alle reads bør komme fra samme område på 16S genet, men altså fra ulike typer bakterier. Ideelt sett skal variasjonen utelukkende skyldes at vi har ulike bakterier i prøven, og ulike bakterier har litt ulik variant av 16S genet. Hele poenget med en slik sekvensering er å få et innblikk i hvor mange typer bakterier vi har i prøven, og hvor ulike de er, samt å gjenkjenne typen basert på at vi vet hvordan dette genet skal se ut for de ulike typene.

La oss sammenligne det vi ser i FastQC for shotgun-data og amplicon-data:

Hvorfor har vi en stor forskjell mellom data-settene når vi ser på Per base sequence content?

Under fanene Sequence Duplication Levels og Overrepresented sequences kan det se ut som at amplicon-dataene er dårlige. Men hva skyldes dette egentlig?



4 R-lab

4.1 Hva er koding?

Koding betyr å lage programmer, det vil si en serie med instruksjoner, som kan tolkes og utføres av en datamaskin. I bunn og grunn kan datamaskiner kun skjønne det som kalles binær-kode. Dette er en serie med 0 og 1 og i praksis helt umulig for oss å lese. Derfor har vi laget programmerings-språk som er slik at de er forståelige for oss mennesker, samtidig som de lar seg oversette eksakt til binær-kode. Dette gjør at vi kan gi datamaskinen instrukser.

Alle programmerings-språk bærer preg av at alle instrukser må være helt eksakt korrekt for at de skal kunne oversettes. Når vi snakker med hverandre kan vi komme unna med ganske upresise utsagn fordi vi har en kontekst som det tolkes inn i, og som gjør at vi skjønner selv ganske upresist språk. Når vi programmerer er det ikke slik. Enhver bitteliten feil gjør at instruksen ikke forstås. I starten blir det mange slike feil. Det er derfor viktig at man er tålmodig og venner seg til at det tar tid å få programmer (noenlunde) feilfrie.

4.2 Google

Dersom du lurer på noe i R, kan det enkleste ofte være å søke på Google. Begynn søket med “R …søkeord…”. Det ligger utrolig mye ute på nettet, og svært ofte vil man finne noe som kan være nyttig. Når man er helt fersk er det selvsagt verre å henge med i forklaringene, men prøv, og lær deg å øke din kompetanse på denne måten. Hvis alle er ferske er det også spesielt effektivt å sitte sammen i grupper, og hyppig utveksle nye oppdagelser av muligheter i R.

4.3 R-minnet

Datamaskiner kan huske. Ekstremt godt. Og det er i grunnen det eneste de kan! Når vi starter RStudio (første gang) så er R-minnet helt tomt. Dersom du tar en titt i Environment fanen i RStudio listes der noe som kalles Global Environment. Her vil det dukke opp noe etterhvert som vi begynner å ta i bruk og opprette objekter i R.

Det vi typisk gjør i ethvert R-program er å opprette en eller flere objekter. Et objekt er simpelthen en bit av R-minnet. Den har alltid et navn, og (nesten) alltid et innhold. La oss nå lage et lite program som leser innholdet i en FASTA-fil inn i R-minnet.

4.4 FASTA-fil

I en øvingsoppgave over ble du bedt om å laste ned fra NCBI-Protein databasen en FASTA-fil med alle sekvenser som matchet et søk etter nitrite reductase (nirK) gener. La oss anta du har lagret dette i fila med navnet nirK_NCBI_protein.fasta (bare bruk ditt navn, dette er kun et eksempel). Det er viktig at du vet eksakt hvor du lagret denne fila på din egen PC.

4.5 Lese inn FASTA-fil

De aller fleste kommandoene vi utfører i R er det vi kaller funksjoner. Vi skal senere se litt på hvordan vi lager slike funksjoner, men først skal vi bare bruke noen ferdige funksjoner. I pakken microseq som vi installerte forrige uke finnes en funksjon med navn readFasta(). Denne kan brukes til å lese en FASTA-fil inn i R-minnet. Her er de 2 linjene med kode som trengs for å få dette til:

library(microseq)
nirK <- readFasta("C:/BIN210/data/NCBI_protein_nirK.fasta")

Ta en titt på filmsnutten les en FASTA-fil med readFasta for å en kort guide som viser hvordan du lager ditt første R-program, og hvordan du kjører programmet.

4.6 Pakker og funksjoner

Alle kommandoer vi bruker i R er funksjoner. En R-pakke er en samling av slike funksjoner. Når du installerer R får du med det vi kaller basis-pakkene, som inneholder alle de grunnleggende funksjonene man alltid trenger. I tillegg kan vi installere ytterligere pakker, noe vi gjorde med blant annet microseq. I koden over brukte vi funksjonen readFasta() for å lese inn en FASTA fil. En slik funksjon finnes ikke i basis R, og vi må da laste inn en R pakke som vi vet har en slik funksjon. Når vi skriver library(microseq) så lastes denne pakken, og alle dens funksjoner blir nå tilgjengelige for oss.

Alle funksjoner i R har en hjelpe-fil, der du kan finne ut mer om hva funksjonen gjør, hvordan den skal brukes og hva den gir tilbake. For å se hjelpe-fila til readFasta() funksjonen så skriver du:

?readFasta

i Console vinduet i RStudio. Da dukker hjelpe-filen opp under Help-fanen. Vi skal mer på slike hjelpe-filer etterhvert.

4.7 Objekter og tilordning

Når vi kjørte program-snutten over fikk vi opprettet et objekt ved navn nirK i R-minnet. Vi bestemmer selv navnet på objekter vi oppretter, men de må bestå av ett ord uten mellomrom. TIPS: Hold deg unna de norske bokstavene æ, ø og å! Selv om det funker noen ganger, så er det slett ikke alltid. Dette gjelder både i navn på objekter, funksjoner og filer. Det er alltid lurt å kalle objekter noe fornuftig som reflekterer det de inneholder.

En tilordning er det vi gjør når vi bruker operatoren <-. I koden over ser vi at vi bruker funksjonen readFasta(), og denne leser inn data fra en fil. Disse dataene tilordner vi til objektet nirK, som opprettes der og da. Operatoren <- likner en pil, og det er et poeng. Den sier at data flyttes fra funksjonen og inn i objektet nirK. Det er faktisk lov å snu denne pila i R, og vi kunne ha skrevet

readFasta("C:/BIN210/data/NCBI_protein_nirK.fasta") -> nirK

noe vi gjør innimellom.

Det er også lov å bruke likhets-tegnet = ved tilordning. Siden dette likner andre språk, velger noen dette. Det anbefales ikke! For det første kan denne operatoren ikke snus, slik pila kan, og dessuten bruker vi likhets-tegnet i en litt annen sammenheng i R. Unngå forvirring, og bruk pila <-.

4.8 Tabeller

R er spesielt godt egnet for å håndtere tabeller av data. Når vi leste inn FASTA filen over, så ble resultatet lagret som en tabell i nirK objektet. Akkurat denne tabellen har to kolonner, en som heter Header og en som heter Sequence. Dette kan vi finne ut på flere måter. En måte er å bare skrive navnet på objektet i Console-vinduet, og return:

nirK
## # A tibble: 9,271 x 2
##    Header                             Sequence                                  
##    <chr>                              <chr>                                     
##  1 ADB12477.1 NirK [Methylomonas sp.~ MDTRHKLNSIDWHLMALIGTLTMALVGCQSQPLVAPTTAEI~
##  2 BAC00912.1 nitrite reductase [Hyp~ MIFHRFGAVVGLMMAALLSVPAAHADQHQQMAQKAAPAADA~
##  3 CAD89521.1 nitrite reductase [Hal~ MLSTTRRRTLQWLGLGGAASLAGCATEAPTTAQSLDGADEP~
##  4 GAP62154.1 nitrite reductase [Ard~ MDYIPDITFTLQTSLSEKGLGFLGVGGDIDGEFNPTLVVPK~
##  5 AKJ30874.1 nitrite reductase [[Po~ MFARVLAPIPFAIAVFVGAAMAGPLTSVAQAAPATAQKEVV~
##  6 CAE46530.1 nitrite reductase [Hal~ MLSTTRRRTLQWLGLGGVASLAGCATKSPTAAQSLDETEEP~
##  7 AYM84367.1 nitrite reductase [Agr~ MTETFHITRRNILAGAAFAGAVAPMVIPSAASAAEEKKAAA~
##  8 AYM65167.1 nitrite reductase [Agr~ MTETFHITRRNILAGAAFAGAIAPMIIPSAAGAAEEKKAAA~
##  9 AYM60637.1 nitrite reductase (pla~ MTETFHITRRNILAGAAFAGAIAPMIIPSAAGAAEEKKTAA~
## 10 AYM08378.1 nitrite reductase [Agr~ MTETFHITRRNILAGAAFAGAVAPMVIPSAASAAEEKKAAA~
## # ... with 9,261 more rows

Da listes hele eller deler av objektet, og vi kan se at dette er en tabell med navn på de to kolonnene.

I RStudio er det også et Environment vindu (nede til venstre?). Her listes alle eksisterende objekter i R minnet. Som du så i filmsnutten over kan vi herfra også åpne tabeller i Viewer-vinduet i RStudio.

I R kalles en tabell enten data.frame eller tibble. Den første er den opprinnelige varianten av en tabell, den andre er nyere, og kom med det som kalles tidyverse-pakkene. Vi installerte også disse sist uke, og skal ta i bruk funksjoner fra tidyverse straks.

Tabellen over er en tibble. Den viser bare fram de første radene, og utskriften under viser at det er mange flere rader. Størrelsen på tabellen kan også avleses i Environment-vinduet, obs. (obervasjoner) er det samme som rader, og variables er kolonner.

4.9 Data-typer

Hva kan en tabell inneholde? I de aller fleste tilfellene vil den enten inneholde tekst eller tall. Det er dette vi kaller data typer. La oss først fokusere på data typen tekst, siden vi skal jobbe mye med sekvenser, og disse oppfattes enklest som tekster.

En tekst i R angis ved at den omsluttes av anførselstegn. Vi kan derfor opprette et objekt og fylle det med en tekst, slik:

tulle_objekt <- "en tekst"

Det er viktig at du skjønner at navnet på objektet her er tulle_objekt (uten anførselstegn) og innholdet er "en tekst" (med anførselstegn). Dette objektet er ikke en tabell, men det vi kaller en vektor. Vi skal se mer på vektorer straks.

Det er flere data typer enn bare tekster (character) og tall (numeric eller integer), og vi skal se disse etterhvert. For mer om data-typer kan du se her. Ikke bare les, men prøv eksemplene selv og eksperimenter!

4.10 Vektorer

En data struktur sier noe om hvordan data er organisert i et objekt. Vi så at nirK ble en tabell, og dette er et eksempel på en slik data struktur. Den mest basale data strukturen er en vektor.

En vektor er en enkel lineær rekke av data, for eksempel noe slikt:

tekst_vektor <- c("en tekst", "enda en tekst", "mer tekst")

Her opprettes objektet tekst_vektor og fylles med en vektor bestående av 3 tekster. Legge merke til hvordan vi bruker funksjonen c() og deretter fyller inn elementene vi vil ha i vektoren, adskilt med komma. Husk at tekster alltid skal omsluttes av anførselstegn.

I en vektor kan vi aldri blande ulike data typer. Det betyr at dersom første elementet er en tekst, så må de andre også være tekster (samme med tall etc.).

En vektor har alltid en lengde, og funksjonen length kan brukes for å se denne

length(tekst_vektor)
## [1] 3

som da er antall elementer i vektoren.

Hver kolonne i en tabell er en vektor. Tabellen nirK har altså to kolonner som begge er vektorer fylt med tekster. I en tabell må alle kolonne-vektorene ha samme lengde. De ulike kolonnene kan imidlertid ha ulike data typer, det vil si at vi nå kan slenge på en ny kolonne i tabellen, og fylle denne med tall. Kravet er at alt i samme kolonne (samme vektor) har samme data type.

For mer om data-strukturer kan du se her: vektorer og tabeller. Her omtales en tabell som en data.frame, men det kan også være en tibble. I de fleste tilfeller spiller det ingen rolle for oss om dataene ligger i en data.frame eller en tibble, tenk på begge som tabeller.

4.11 Manipulere tabeller

La oss fortsette på scriptet vi startet over, der vi leste inn en FASTA fil. Det kunne vært fint å fått beregnet hvor lange de ulike sekvensene er, altså hvor mange tegn det er i hvert protein. Da vil det være naturlig at vi oppbevarer dette som en ny kolonne i den samme tabellen. Her er en linje du kan legge til scriptet:

nirK <- mutate(nirK, Length = str_length(Sequence))

Dersom du nå kjører hele scriptet på nytt (bruk Source-knappen i RStudio), så skal tabellen nirK ha fått en ny kolonne med navn Length.

Her brukte vi funksjonen mutate(). Denne tar som input an tabell (nirK) og muterer denne, dvs

  • enten legger til en ny kolonne
  • eller endrer innholdet i en eksistrende kolonne

Vi sier til funksjonen at Length = str_length(Sequence). Til venstre for likhetstegnet angir vi et kolonne-navn. Siden dette navnet ikke finnes i tabellen fra før, betyr det at vi vil ha en ny kolonne med dette navnet. Hadde vi angitt en eksisterende kolonne, så hadde den nå fått nytt innhold. Innholdet er uansett det som genereres på høyresiden av likhetstegnet. Her bruker vi funksjonen str_length(), som teller antall tegn i tekster. Vi sender hele kolonnen Sequence som input til denne, og får talt opp antall tegn i hver av tekstene. Funksjonen str_length() må altså få en vektor av tekster som input, og leverer tilbake en vektor med heltall (antall tegn). Vi kan godt prøve den på vår tekst_vektor fra tidligere (antar den fortsatt ligger i R minnet):

str_length(tekst_vektor)
## [1]  8 13  9

Input trenger altså ikke være en kolonne i en tabell, enhver vektor fungerer fint.

Merk at vi her brukte likhetstegnet, ikke pila <-. Generelt kan vi si at når vi angir tilordninger i en input til en funksjon (her mutate()), så bruker vi alltid =. Faktisk vil bruk av <- her gi en litt annen oppførsel (prøv).

Merk også at tabellen nirK opptrer på begge sider av tilordningen <-. På høyre siden sender vi tabellen som input til mutate(), mens vi også sier at resultatet skal tilordnes nirK på venstre siden. Dette skal leses omtrent slik: Send eksisterende versjon av nirK som input til mutate(), og når denne funksjonen er ferdig, ta resultatet og legg det som nytt (oppdatert) innhold i objektet nirK. Det er derfor ganske vanlig å se samme objekt på begge sider av en tilordning, det betyr at objektets innhold skal oppdateres, delvis basert på det gamle innholdet.

4.12 Operatoren %>%

La oss introdusere rørgate operatoren %>%. Den kalles pipe på engelsk, og er vanlig i R programmering. Vi bruker denne når vi vil serie-koble flere operasjoner. Ofte vil vi utføre en serie av operasjoner på et objekt. Dette kan gjøres ved at vi sender objektet som input til en funksjon, og får tilbake en oppdatert versjon, som lagres. Så sendes denne som input til neste funksjon, og vi får en ny versjon, som lagres.

Med rørgate operatoren slipper vi mellomlagringen, og koden blir også enklere å lese, som en ‘strøm’ av operasjoner. Her er scriptet i en ny versjon, der jeg har tatt i bruk rørgate:

library(microseq)
readFasta("C:/BIN210/data/NCBI_protein_nirK.fasta") %>% 
  mutate(Length = str_length(Sequence)) -> nirK

Som før leser vi inn FASTA filen, men resultatet lagres ikke i noe objekt, det sendes direkte gjennom rørgaten til neste funksjon, mutate(). Legg også merke til at første input til mutate() nå ikke spesifiseres (sammenlign med koden over). Det er fordi denne nå kommer gjennom rørgaten. Det er alltid første input til funksjonen som kommer gjennom rørgaten. Til slutt tilordner vi til nirK objektet ved å bruke pila snudd mot høyre.

Både den opprinnelige og denne formen for koding vil fungere, scriptet gjør eksakt det samme. Men, den siste måten er oftest hakket raskere og dessuten enklere å lese(?). Vi ser her at det er en strøm av operasjoner (vel, bare to) som til slutt ender med objektet nirK. Vi vil se mye bruk av dette framover.

På tastaturet får du fram rørgate operatoren med SHIFT-CTRL-M.

4.13 Histogram over sekvens-lengder

La oss avslutte med å lage en figur.

Det er mange funksjoner i R for å lage figurer. Vi vil imidlertid bruke funksjoner som gjør bruk at pakken ggplot2, som er en del at tidyverse. Funksjonen ggplot() er sentral. Denne lager et grafisk objekt, som vi kan legge til egenskaper hos. Bokstavene GG står for Grammar of Graphics. Måten ggplot() brukes på gjennomsyres av denne filosofien. I BIN210 vil vi bare ta dette i bruk, og så vil forhåpentligvis ting klarne litt underveis. Det skal sies at ggplot-universet er stort, og det vil ta noe tid før man har dette under huden. Det har imidlertid blitt en slags standard for plotting av vitenskapelige figurer, og vel verdt litt innsats på sikt.

Her er scriptet etter at vi har lagt til kode for å lage et histogram over sekvens-lengdene:

library(tidyverse)
library(microseq)

readFasta("C:/BIN210/data/NCBI_protein_nirK.fasta") %>% 
  mutate(Length = str_length(Sequence)) -> nirK

fig <- ggplot(nirK) +
  geom_histogram(aes(x = Length), bins = 50)
print(fig)

Først, vi må nå også laste inn tidyverse for at ggplot() skal bli tilgjengelig for oss.

Legg merke til at vi sender tabellen nirK som input til ggplot(), det er (nesten) alltid en tabell som er input til denne funksjonen. Deretter kommer en +. Her betyr det ‘legg til et nytt lag på figuren’. Det minner litt om rørgaten vi så over, men operatoren + har denne spesielle betydningen kun i forbindelse med ggplot() (den betyr ellers vanlig pluss, av tall).

Det nye laget vi legger på figuren er det vi kaller et geom. Det finnes drøssevis av slike geomer i ggplot2-pakken, og geom_histogram() lager altså et histogram. Alle slike geomer må få angitt en estetikk, med funksjonen aes() (engelsk aesthetics). Her angis hvilke kolonner fra tabellen som skal brukes. Et histogram må ha angitt hva vi skal ha på x-aksen, og vi vil bruke kolonnen Length. Deretter angir vi bins = 50. Et histogram deler alltid x-aksen opp i intervaller, eller bins på engelsk. Vi vil altså ha delt dette opp i 50 slike intervaller. Dette påvirker stolpe-bredden i histogrammet (jo flere bins jo flere og smalere stolper).

Fra histogrammet ser vi at de aller fleste nirK proteinene er litt under 400 tegn lange, men at det også er en ny topp rundt 500. Det kan muligens bety at vi egentlig har to litt ulike proteiner her, selv om begge ble fanget opp av søket vårt hos NCBI? Slike tekst-søk vil alltid være litt upresise, så dette er ikke så overraskende egentlig.

4.14 På nettet

Det finnes enorme mengder med informasjon på nettet omkring bruk av R, og tidyverse pakkene. Dersom du vil lære mer kan du for eksempel ta en titt på DataCamp kurs på nettet. De har noen gratis introduksjoner, som ‘lokkemat’. Vi anbefaler at dere tar en titt på disse.

Det finnes ellers enormt mye gratis på nettet, del med alle på Canvas dersom du finner noe du vil anbefale.