Informatie

Wat is de lengte van het gen bij het berekenen van TPM (transcripten per miljoen)


Wat is de lengte van het gen bij het berekenen van TPM (transcripten per miljoen)? Stel dat ik een datasetmatrix heb met k rijen (elke rij is een gen) en n kolommen (elke kolom is een voorbeeld), is er een manier om de dataset over te zetten naar TPM? Erg bedankt!


Stel dat we hebben

$$ : : egin{bmatrix} Gene : Name & Rep1 : Count & Rep2 : Count & Rep3 : Count A : (2kb) & 10 & 15 & 30 B: ( 5kb) & 25 & 25 & 70 C: (2kb) & 6 & 10 & 17 D: (8kb) & 1 & 1 & 3 end{bmatrix} : :$$


Om dit om te zetten in TPM-indeling, we moeten normaliseren voor genlengte, en normaliseer vervolgens voor gendiepte, in die volgorde.


Normaliseren voor genlengte...

  • Stap(pen) om uit te voeren: deel elke herhalingstelling door de lengte van het respectieve gen.

Voor $Gene A$ heeft het een lengte van $2kb$ (kilobasen), met replica's van $10, 15,$ en $30$. Als u deze bewerking op de hele tafel uitvoert,

$$ : : egin{bmatrix} Gene : Naam & Rep1 : Count & Rep2 : Count & Rep3 : Count A : (2kb) & frac{10}{2} & frac {15}{2} & frac{30}{2} B: (5kb) & frac{25}{5} & frac{25}{5} & frac{70}{5} C: (2kb) & frac{6}{2} & frac{10}{2} & frac{17}{2} D: (8kb) & frac{1}{ 8} & frac{1}{8} & frac{3}{8} end{bmatrix} : :$$

opbrengsten

$$ : : egin{bmatrix} Gene : Naam & Rep1 : Count & Rep2 : Count & Rep3 : Count A : (2kb) & 5 & 7.5 & 15 B: ( 5kb) & 5 & 5 & 35 C: (2kb) & 3 & 5 & 8.5 D: (8kb) & .125 & .125 & .375 end{bmatrix} : :$$


Normaliseren voor gendiepte...

  • Stap(pen) die u moet uitvoeren: 1) Tel alle tellingen op binnen elke herhalingskolom; 2) Deel elke kolomsom door de gewenste diepte (dit levert schaalfactoren op); 3) Deel elke herhalingstelling binnen een kolom door de respectieve schaalfactor.

$Rep1$

  • Som: $5 + 5 + 3 + .125 = 13.125$
  • Schaalfactor: $frac{13.125}{1.000.000} = 1.325 imes 10^{-5}$

$Rep2$

  • Som: $7,5 + 5 + 5 + 0,125 = 17,625 $
  • Schaalfactor: $frac{17.625}{1.000.000} = 1.7625 imes 10^{-5}$

$Rep3$

  • Som: $15 + 35 + 8.5 + .375 = 58.875$
  • Schaalfactor: $frac{13.125}{1.000.000} = 5.8875 imes 10^{-5}$

Door deze waarden op de tabel toe te passen,

$$ : : egin{bmatrix} Gene : Name & Rep1 : Count & Rep2 : Count & Rep3 : Count A : (2kb) & frac{5}{1.325 imes 10^ {-5}} & frac{7.5}{1.7625 imes 10^{-5}} & frac{15}{5.8875 imes 10^{-5}} B: (5kb) & frac {5}{1.325 imes 10^{-5}} & frac{5}{1.7625 imes 10^{-5}} & frac{35}{5.8875 imes 10^{-5}} C: (2kb) & frac{3}{1.325 imes 10^{-5}} & frac{5}{1.7625 imes 10^{-5}} & frac{8.5}{5.8875 imes 10^{-5}} D: (8kb) & frac{.125}{1.325 imes 10^{-5}} & frac{.125}{1.7625 imes 10^{-5} } & frac{.375}{5.8875 imes 10^{-5}} end{bmatrix} : : $$

opbrengsten

$$ : : egin{bmatrix} Gene : Naam & Rep1 : Count & Rep2 : Count & Rep3 : Count A : (2kb) & 377358.49 & 425531.91 & 254777.07 B: ( 5kb) & 377358.49 & 283687.94 & 594479.83 C: (2kb) & 226415.09 & 283687.94 & 144373.67 D: (8kb) & 9433.96 & 7092.20 & 6369.43 end{bmatrix} : :$$

De dataset is nu TPM-geformatteerd, wat een eenvoudigere analyse van afgelezen verhoudingen in het hele monster biedt.


RNA-Seq-normalisatie uitgelegd

RNA-Seq (afkorting van RNA-sequencing) is een soort experiment waarmee we genexpressie kunnen meten. De sequencingstap produceert een groot aantal (tientallen miljoenen) cDNA 1-fragmentsequenties die reads worden genoemd. Elke aflezing vertegenwoordigt een deel van een RNA-molecuul in het monster 2 .

Vervolgens wijzen we elke read toe aan een van de isovormen en tellen we hoeveel reads elke isovorm heeft.

Als al het andere gelijk is, geldt dat hoe overvloediger een isovorm is, hoe meer fragmenten ervan waarschijnlijk worden gesequenced. Daarom kunnen we leestellingen gebruiken als een proxy voor isovorm-abundantie.

Omdat "al het andere" echter nooit gelijk is, moeten de tellingen worden aangepast om vergelijkbaar te zijn tussen isovormen, monsters en experimenten. Hier zullen we deze aanpassingen onderzoeken en waarom ze nodig zijn.

Overweeg de volgende in kaart gebrachte leest van een RNA-Seq-experiment. Welke isovorm komt meer voor, de rode of de gele?

De gele isovorm heeft meer aflezingen, maar is ook veel langer dan de rode. Hoe langer de isovorm, hoe meer fragmenten (en dus leest) we zouden verwachten dat deze zal genereren.

Om leestellingen over isovormen te kunnen vergelijken, delen we de tellingen door de isovormlengte. Het is ook gebruikelijk om het getal te vermenigvuldigen met (1000) , waardoor leest per kilobase:

waarbij (n_i) het aantal reads is dat is toegewezen aan isovorm (i) , en (l_i) de lengte van die isovorm is.


Genexpressie-eenheden uitgelegd: RPM, RPKM, FPKM en TPM

Bij de analyse van RNA-seq-genexpressiegegevens komen we verschillende expressie-eenheden tegen, zoals RPM, RPKM, FPKM en onbewerkte reads. Meestal is het moeilijk om de onderliggende basismethodologie te begrijpen om deze eenheden te berekenen uit in kaart gebrachte sequentiegegevens. Hier heb ik geprobeerd deze eenheden op een veel eenvoudigere manier uit te leggen.

Waarom verschillende genormaliseerde expressie-eenheden:

De expressie-eenheden bieden een digitale maatstaf voor de overvloed aan transcripten. Genormaliseerde expressie-eenheden zijn nodig om technische vooroordelen in sequentiegegevens te verwijderen, zoals de diepte van de sequencing (meer sequencing-diepte produceert meer leestellingen voor genen die op hetzelfde niveau tot expressie worden gebracht) en genlengte (verschillen in genlengte genereren ongelijke leestellingen voor genen die op hetzelfde niveau tot expressie worden gebracht) niveau langer het gen meer de gelezen telling).

Genexpressie-eenheden en berekening:

RPM (Reads per miljoen toegewezen reads)

U hebt bijvoorbeeld één bibliotheek gesequenced met 5 miljoen (M) leesbewerkingen. Onder hen kwam in totaal 4 M overeen met de genoomsequentie en 5000 reads die overeenkwamen met een bepaald gen.

  • RPM houdt geen rekening met de normalisatie van de transcriptlengte.
  • RPM Geschikt voor sequencing-protocollen waarbij uitlezingen worden gegenereerd ongeacht de genlengte

RPKM (Leest per kilo basis per miljoen toegewezen reads)

Hier normaliseert 103 voor genlengte en 106 voor sequentiedieptefactor.

FPKM (Fragmenten per kilobasis per miljoen toegewezen reads) is analoog aan RPKM en wordt vooral gebruikt in gepaarde RNA-seq-experimenten. In gepaarde RNA-seq-experimenten worden twee (links en rechts) aflezingen gesequenced van hetzelfde DNA-fragment. Wanneer we gepaarde-end-gegevens in kaart brengen, kunnen beide reads of slechts één read met hoge kwaliteit van een fragment worden toegewezen aan de referentiesequentie. Om verwarring of meervoudig tellen te voorkomen, worden de fragmenten waaraan beide of enkelvoudige leestoewijzing zijn toegewezen, geteld en weergegeven voor FPKM-berekening.

U hebt bijvoorbeeld één bibliotheek gesequenced met 5 miljoen leesbewerkingen. Onder hen kwam in totaal 4 M overeen met de genoomsequentie en 5000 reads die overeenkwamen met een bepaald gen met een lengte van 2000 bp.

  • RPKM houdt rekening met de genlengte voor normalisatie
  • RPKM is geschikt voor sequencing-protocollen waarbij reads-sequencing afhangt van de genlengte
  • Gebruikt in single-end RNA-seq-experimenten (FPKM voor gepaarde RNA-seq-gegevens)

TPM (transcriptie per miljoen)

Hier verwijst leeslengte naar het gemiddelde aantal nucleotiden dat aan een gen is toegewezen.


Wat is de lengte van het gen bij het berekenen van TPM (transcripten per miljoen) - Biologie

Bij de analyse van RNA-seq-genexpressiegegevens komen we verschillende expressie-eenheden tegen, zoals RPM, RPKM, FPKM en onbewerkte uitlezingen. Meestal is het moeilijk om de onderliggende basismethodologie te begrijpen om deze eenheden te berekenen uit in kaart gebrachte sequentiegegevens.

Ik heb veel berichten gezien over dergelijke normalisatievragen en hun verwarring onder lezers. Daarom heb ik hier geprobeerd deze eenheden op een veel eenvoudigere manier uit te leggen (waarbij ik complexe wiskundige uitdrukkingen heb vermeden).

Waarom verschillende genormaliseerde expressie-eenheden:

De expressie-eenheden bieden een digitale maatstaf voor de overvloed aan transcripten. Genormaliseerde expressie-eenheden zijn nodig om technische vooroordelen in sequentiegegevens te verwijderen, zoals de diepte van de sequencing (meer sequencing-diepte produceert meer leestellingen voor genen die op hetzelfde niveau tot expressie worden gebracht) en genlengte (verschillen in genlengte genereren ongelijke leestellingen voor genen die op hetzelfde niveau tot expressie worden gebracht) niveau langer het gen meer de gelezen telling).

Genexpressie-eenheden en berekening:

1. RPM (Leesbewerkingen per miljoen toegewezen leesbewerkingen)

U hebt bijvoorbeeld één bibliotheek gesequenced met 5 miljoen (M) leesbewerkingen. Onder hen kwam in totaal 4 M overeen met de genoomsequentie en 5000 reads die overeenkwamen met een bepaald gen.

  • RPM houdt geen rekening met de normalisatie van de transcriptlengte.
  • RPM Geschikt voor sequencing-protocollen waarbij uitlezingen worden gegenereerd, ongeacht de genlengte! [voer afbeeldingsbeschrijving in

2. RPKM (leeswaarden per kilobasis per miljoen toegewezen aflezingen)

Hier normaliseert 10 ^ 3 voor genlengte en 10 ^ 6 voor sequencing-dieptefactor.

FPKM (Fragments per kilo base per miljoen toegewezen reads) is analoog aan RPKM en wordt vooral gebruikt in gepaarde RNA-seq-experimenten. In gepaarde RNA-seq-experimenten worden twee (links en rechts) aflezingen gesequenced van hetzelfde DNA-fragment. Wanneer we gepaarde-end-gegevens in kaart brengen, kunnen beide reads of slechts één read met hoge kwaliteit van een fragment worden toegewezen aan de referentiesequentie. Om verwarring of meervoudig tellen te voorkomen, worden de fragmenten waaraan beide of enkelvoudige leestoewijzing zijn toegewezen, geteld en weergegeven voor FPKM-berekening.

U hebt bijvoorbeeld één bibliotheek gesequenced met 5 miljoen leesbewerkingen. Onder hen kwam in totaal 4 M overeen met de genoomsequentie en 5000 reads die overeenkwamen met een bepaald gen met een lengte van 2000 bp.


RNA-Seq-gegevens schalen en normaliseren

Transcriptlengtes

Reads (Fragments) Per Kilobase Million (RPKM) en Transcripts Per Million (TPM) zijn statistieken om genexpressie te schalen om twee doelen te bereiken

  1. Maak de expressie van genen vergelijkbaar tussen monsters.
  2. Maak de expressie van verschillende genen vergelijkbaar.

Voor (1.) zullen de bibliotheekgroottes (aantal totale leesbewerkingen) altijd verschillen tussen monsters als een technisch artefact van RNA-sequencing. Voor (2.) is de grootte van RNA-transcripten van elk gen anders en we verwachten dat er meer reads worden geteld in grotere transcripten.

RPKM en TPM zijn zeer vergelijkbare statistieken. Voor elk gen in elk monster...

[ TPM_ = frac><>/1000>>^G frac<>><>/1000>> imes 1E6 ] Waarbij (G) het totale aantal genen is. Het verschil is subtiel, maar merk op dat de bibliotheekgrootte voor RPKM we eerst de bibliotheekgrootte schalen (som van onbewerkte tellingen), waar we voor TPM eerst schalen voor de transcriptgrootte en vervolgens schalen met de som van deze getransformeerde tellingen. Alleen TPM zorgt ervoor dat de geschaalde bibliotheekgroottes gelijk zijn voor alle samples, waarbij de som van de RPKM-waarden tussen samples verschilt.


Wat is de lengte van het gen bij het berekenen van TPM (transcripten per miljoen) - Biologie

Loop string stropdas vanaf de opdrachtregel als volgt:
De belangrijkste invoer van het programma is een BAM-bestand met RNA-Seq-leestoewijzingen die moeten worden gesorteerd op hun genomische locatie (bijvoorbeeld de geaccepteerd_hits.bam bestand geproduceerd door TopHat of de uitvoer van HISAT2 na het sorteren en converteren met samtools zoals hieronder uitgelegd).

Invoerbestanden

Het bestand is het resultaat van de bovenstaande opdracht (alns.sorted.bam) kan worden gebruikt als invoer voor StringTie.

Elke gesplitste leesuitlijning (d.w.z. een uitlijning over ten minste één knooppunt) in het SAM-invoerbestand moet de tag bevatten XS om de genomische streng aan te geven die het RNA produceerde waarvan de sequentie werd bepaald. Uitlijningen geproduceerd door TopHat en HISAT2 (wanneer uitgevoerd met de optie --dta) bevatten deze tag al, maar als u een andere leesmapper gebruikt, moet u controleren of deze XS-tag is opgenomen voor gesplitste uitlijningen.

OPMERKING: zorg ervoor dat u HISAT2 uitvoert met de optie --dta voor uitlijning, anders zullen uw resultaten eronder lijden.

Als optie kan een referentie-annotatiebestand in GTF/GFF3-formaat aan StringTie worden verstrekt. In dit geval zal StringTie er de voorkeur aan geven deze "bekende" genen uit het annotatiebestand te gebruiken, en voor degenen die worden uitgedrukt, berekent het de dekking, TPM en FPKM-waarden. Het zal ook extra transcripties produceren om rekening te houden met RNA-seq-gegevens die niet worden gedekt door (of verklaard door) de annotatie. Merk op dat als optie -e niet wordt gebruikt, de referentietranscripties volledig moeten worden gedekt door reads om te worden opgenomen in de uitvoer van StringTie. In dat geval worden ook andere transcripties afgedrukt die door StringTie zijn samengesteld uit de gegevens en niet aanwezig zijn in het referentiebestand.

OPMERKING: we raden u ten zeerste aan om annotaties op te geven als u een genoom analyseert dat goed is geannoteerd, zoals mensen, muizen of andere modelorganismen.

Uitvoerbestanden

  1. De belangrijkste output van Stringtie is een GTF-bestand met de verzamelde transcripties
  2. Gen-abundanties in door tabs gescheiden formaat
  3. Volledig gedekte transcripties die overeenkomen met de referentieannotatie, in GTF-indeling
  4. Bestanden (tabellen) die nodig zijn als invoer voor Ballgown, die ze gebruikt om differentiële expressie te schatten
  5. In de samenvoegmodus, een samengevoegd GTF-bestand uit een set GTF-bestanden

Beschrijving van de waarden van elke kolom:

  • volgnaam: Geeft het chromosoom, contig of scaffold voor dit transcript aan. Hier bevindt het geassembleerde transcript zich op chromosoom X.
  • bron: De bron van het GTF-bestand. Aangezien dit voorbeeld is gemaakt door StringTie, toont deze kolom gewoon 'StringTie'.
  • functie: Kenmerktype, bijv. exon, transcript, mRNA, 5'UTR).
  • begin: Startpositie van de functie (exon, transcript, etc), met behulp van een 1-gebaseerde index.
  • einde: Eindpositie van het kenmerk, met behulp van een 1-gebaseerde index.
  • scoren: Een betrouwbaarheidsscore voor het samengestelde transcript. Momenteel wordt dit veld niet gebruikt, en StringTie rapporteert een constante waarde van 1000 als het transcript een verbinding heeft met een leesuitlijningsbundel.
  • strand: Als het transcript zich op de voorwaartse streng bevindt, '+'. Als het transcript zich op de omgekeerde streng bevindt, '-'.
  • kader: Frame of fase van CDS-functies. StringTie gebruikt dit veld niet en registreert gewoon een ".".
  • attributen: Een door puntkomma's gescheiden lijst van tag-waardeparen, met aanvullende informatie over elk kenmerk. Afhankelijk van of een instantie een transcript of een exon is en of het transcript overeenkomt met het referentieannotatiebestand dat door de gebruiker is verstrekt, zal de inhoud van het attributenveld verschillen. De volgende lijst beschrijft de mogelijke attributen die in deze kolom worden weergegeven:
    • gene_id: Een unieke identifier voor een enkel gen en zijn onderliggende transcript en exons op basis van de bestandsnaam van de uitlijningen.
    • transcript_id: Een unieke identifier voor een enkel transcript en de onderliggende exons op basis van de bestandsnaam van de uitlijningen.
    • exon_number: Een unieke identifier voor een enkele exon, beginnend bij 1, binnen een bepaald transcript.
    • reference_id: de transcript_id in de referentie-annotatie (optioneel) waarmee de instantie overeenkomt.
    • ref_gene_id: de gene_id in de referentieannotatie (optioneel) waarmee de instantie overeenkwam.
    • ref_gene_name: de gene_name in de referentieannotatie (optioneel) waarmee de instantie overeenkwam.
    • cov: de gemiddelde dekking per base voor het transcript of exon.
    • FPKM: Fragmenten per kilobase transcript per miljoen gelezen paren. Dit is het aantal leesparen dat overeenkomt met deze functie, genormaliseerd door het totale aantal fragmenten waarvan de sequentie is bepaald (in miljoenen) en de lengte van het transcript (in kilobasen).
    • TPM: transcripties per miljoen. Dit is het aantal transcripten van dit specifieke gen, eerst genormaliseerd op genlengte en vervolgens op sequentiediepte (in miljoenen) in het monster. Een gedetailleerde uitleg en een vergelijking van TPM en FPKM vindt u hier, en TPM is hier gedefinieerd door B. Li en C. Dewey.
    • Kolom 1 / Gen-ID: De genidentificatie komt van de referentieannotatie die bij de -G-optie wordt geleverd. Als er geen verwijzing is opgegeven, wordt dit veld vervangen door het naamvoorvoegsel voor uitvoertranscripten (-l).
    • Kolom 2 / Gennaam: Dit veld bevat de gennaam in de referentieannotatie die bij de -G-optie wordt geleverd. Als er geen referentie is opgegeven, wordt dit veld gevuld met '-'.
    • Kolom 3 / Verwijzing: Naam van de referentiesequentie die werd gebruikt bij de uitlijning van de uitlezingen. Gelijk aan de 3e kolom in de .SAM-uitlijning.
    • Kolom 4 / Strand: '+' geeft aan dat het gen op de voorwaartse streng zit, '-' voor de achterwaartse streng.
    • Kolom 5 / Begin: Startpositie van het gen (1-gebaseerde index).
    • Kolom 6 / Einde: Eindpositie van het gen (1-gebaseerde index).
    • Kolom 7 / Dekking: Per-base dekking van het gen.
    • Kolom 8 / FPKM: genormaliseerd expressieniveau in FPKM-eenheden (zie vorige paragraaf).
    • Kolom 9 / TPM: genormaliseerd expressieniveau in RPM-eenheden (zie vorige paragraaf).

    3. Volledig gedekte transcripties die overeenkomen met de transcripten van de referentieannotatie (in GTF-formaat)

    Als StringTie wordt uitgevoerd met de optie -C <cov_refs.gtf> (vereist -G <reference_annotation> ), retourneert het een bestand met alle transcripten in de referentieannotatie die volledig zijn gedekt, van begin tot eind, door reads. Het uitvoerformaat is een GTF-bestand zoals hierboven beschreven. Elke regel van de GTF komt overeen met een gen of transcript in de referentieannotatie.

    4. Ballgown-invoertabelbestanden

    Als StringTie wordt uitgevoerd met de optie -B, retourneert het een Ballgown-invoertabelbestand, dat dekkingsgegevens voor alle transcripties bevat. De uitvoertabelbestanden worden in dezelfde map geplaatst als de hoofd-GTF-uitvoer. Deze tabellen hebben de volgende specifieke namen: (1) e2t.ctab, (2) e_data.ctab, (3) i2t.ctab, (4) i_data.ctab en (5) t_data.ctab. Een gedetailleerde beschrijving van elk van deze vijf vereiste inputs voor Ballgown is te vinden op de Ballgown-site, via deze link.

    5. Samenvoegmodus: samengevoegde GTF

    Als StringTie wordt uitgevoerd met de optie --merge, neemt het als invoer een lijst van GTF/GFF-bestanden en voegt het deze transcripten samen/assembleert het tot een niet-redundante set transcripties. Deze stap creëert een uniforme set transcripten voor alle monsters om de stroomafwaartse berekening van differentieel uitgedrukte niveaus voor alle transcripten onder de verschillende experimentele omstandigheden te vergemakkelijken. De uitvoer is een samengevoegd GTF-bestand met alle samengevoegde genmodellen, maar zonder numerieke resultaten op dekking, FPKM en TPM. Met deze samengevoegde GTF kan StringTie de abundanties opnieuw schatten door het opnieuw uit te voeren met de optie -e op de originele set uitlijningsbestanden, zoals geïllustreerd in de onderstaande afbeelding.

    Transcriptsamenstellingen evalueren

    Een eenvoudige manier om meer informatie te krijgen over de transcripten die door StringTie zijn verzameld (samenvatting van het aantal genen en transcripten, nieuw versus bekend, enz.), of zelfs basistracking van geassembleerde isovormen over meerdere RNA-Seq-experimenten uit te voeren, is door het gffcompare-programma te gebruiken . Basisgebruiksinformatie en downloadopties voor dit programma zijn te vinden op de GFF-hulpprogramma's-pagina.

    Differentiële expressieanalyse

    1. wijs voor elk RNA-Seq-monster de uitlezingen toe aan het genoom met HISAT2 met behulp van de --dta optie. Het wordt ten zeerste aanbevolen om de referentie-annotatie-informatie te gebruiken bij het in kaart brengen van de uitlezingen, die ofwel kunnen worden ingebed in de genoomindex (gebouwd met de --ss en --exon opties, zie HISAT2-handleiding), of apart geleverd tijdens runtime (met behulp van de --bekend-splicesite-infile optie van HISAT2). De SAM-uitvoer van elke HISAT2-run moet worden gesorteerd en geconverteerd naar BAM met behulp van samtools zoals hierboven uitgelegd.
    2. Voer voor elk RNA-Seq-monster StringTie uit om de leesuitlijningen te verzamelen die in de vorige stap zijn verkregen. Het wordt aanbevolen om StringTie uit te voeren met de optie -G als de referentieannotatie beschikbaar is.
    3. voer StringTie uit met --samenvoegen om een ​​niet-redundante set transcripten te genereren die zijn waargenomen in een van de eerder geassembleerde RNA-Seq-monsters. De stringtie --merge mode neemt als invoer een lijst van alle verzamelde transcriptiebestanden (in GTF-formaat) die eerder voor elk monster zijn verkregen, evenals een referentie-annotatiebestand (-G optie) indien beschikbaar.
    4. voer voor elk RNA-Seq-monster StringTie uit met behulp van de -B/-b en -e opties om transcript-abundanties te schatten en leesdekkingstabellen voor Ballgown te genereren. De optie -e is niet vereist, maar wordt aanbevolen voor deze run om nauwkeurigere schattingen van de hoeveelheid van de invoertranscripten te produceren. Elke StringTie-run in deze stap zal de gesorteerde leesuitlijningen (BAM-bestand) die in stap 1 zijn verkregen voor het overeenkomstige voorbeeld en de -G optie met de samengevoegde transcripties (GTF-bestand) gegenereerd door stringtie --merge in stap 3. Houd er rekening mee dat dit het enige geval is waarbij de -G optie wordt niet gebruikt met een referentieannotatie, maar met de globale, samengevoegde set transcripties zoals waargenomen in alle monsters. (Deze stap is het equivalent van de Tafelmaker stap beschreven in de originele Ballgown-pijplijn.)
    5. Ballgown kan nu worden gebruikt om de in de vorige stap gegenereerde dekkingstabellen te laden en verschillende statistische analyses uit te voeren voor differentiële expressie, plots te genereren, enz.

    Een alternatieve, snellere workflow voor differentiële expressie-analyse kan worden nagestreefd als er geen interesse is in nieuwe isovormen (dwz geassembleerde transcripten die aanwezig zijn in de monsters maar ontbreken in de referentieannotatie), of als alleen een bekende reeks van belang zijnde transcripten het doelwit zijn van de analyse. Dit vereenvoudigde protocol heeft slechts 3 stappen (hieronder afgebeeld) omdat het de individuele assemblage van elk RNA-Seq-monster en de "transcriptie samenvoegen" stap. Deze vereenvoudigde workflow probeert de expressie van een bekende set transcripten direct te schatten en te analyseren, zoals gegeven in de referentie annotatie het dossier.

    Het R-pakket IsovormSwitchAnalyserenR kan worden gebruikt om gennamen toe te kennen aan transcripten die door StringTie zijn samengesteld, wat vooral handig kan zijn in gevallen waarin StringTie deze opdracht niet eenduidig ​​kon uitvoeren.
    De functie importIsoformExpression() + importRdata() van het pakket kan worden gebruikt om de expressie- en annotatiegegevens in R te importeren. Tijdens deze import zal het pakket proberen om waar mogelijk isoform-annotaties op te schonen en te herstellen. Van het resulterende switchAnalyzeRlist-object, IsovormSwitchAnalyserenR kan isoform-switches detecteren samen met voorspelde functionele gevolgen. De functie extractGeneExpression() kan worden gebruikt om een ​​matrix voor genexpressie (lees telling) op te halen voor analyse met andere tools.
    Meer informatie en codevoorbeelden vindt u hier.

    StringTie gebruiken met DESeq2 en edgeR

    DESeq2 en edgeR zijn twee populaire Bioconductor-pakketten voor het analyseren van differentiële expressie, die als invoer een matrix van leestellingen nemen die zijn toegewezen aan bepaalde genomische kenmerken (bijv. Genen). We bieden een Python-script (prepDE.py, of de Python 3-versie: prepDE.py3 ) die kan worden gebruikt om deze gelezen telling-informatie rechtstreeks uit de bestanden te extraheren die zijn gegenereerd door StringTie (uitvoeren met de -e parameter).

    prepDE.py leidt hypothetische leestellingen voor elk transcript af van de dekkingswaarden geschat door StringTie voor elk transcript, met behulp van deze eenvoudige formule: reads_per_transcript = dekking * transcript_len / read_len

    • een optie is om een ​​pad op te geven naar een map die alle voorbeeldsubmappen bevat, met dezelfde structuur als de baljurk map in het StringTie-protocoldocument ter voorbereiding op Ballgown. Standaard (geen -l optie), gaat het script in de huidige map zoeken naar alle submappen met .gtf-bestanden, zoals in dit voorbeeld:
    • Als alternatief kan men een tekstbestand opgeven met voorbeeld-ID's en hun respectieve paden (voorbeeld_lst.txt).

    Gebruik: prepDE.py [opties]
    genereert twee CSV-bestanden met de telmatrices voor genen en transcripten, met behulp van de dekkingswaarden die worden gevonden in de uitvoer van stringtie -e

    Opties:
    -h, --help toon dit helpbericht en sluit af
    -i INGANG, --input=INPUT, --in=INPUTeen map met alle voorbeeld-submappen, of een tekstbestand met voorbeeld-ID en pad naar het GTF-bestand op elke regel [standaard: . ]
    -g Gwaar de genentellingsmatrix moet worden uitgevoerd [standaard: gene_count_matrix.csv]
    -t Twaar de transcript-tellingsmatrix moet worden uitgevoerd [standaard: transcript_count_matrix.csv]
    -l LENGTE, --lengte=LENGTEde gemiddelde leeslengte [standaard: 75]
    -p PATROON, --patroon=PATROONeen reguliere expressie die de voorbeeldsubmappen selecteert
    -c, --clusterof genen die overlappen met verschillende gen-ID's moeten worden geclusterd, waarbij genen met een gen-ID-patroon worden genegeerd (zie hieronder)
    -s STRING, --string=STRINGals een ander voorvoegsel wordt gebruikt voor geneID's die zijn toegewezen door StringTie [standaard: MSTRG]
    -k SLEUTEL, --key=KEYbij clustering, welk voorvoegsel moet worden gebruikt voor gen-ID's die door dit script zijn toegewezen [standaard: prepG]
    --legend=LEGENDEals clustering, waar de transcripties van de legendabestandstoewijzing moeten worden uitgevoerd naar toegewezen gen-ID's [standaard: legend.csv]

    Deze telmatrices (CSV-bestanden) kunnen vervolgens in R worden geïmporteerd voor gebruik door DESeq2 en edgeR (met behulp van de DESeqDataSetFromMatrix en DGELlijst functies, respectievelijk).

    Protocol: StringTie gebruiken met DESeq2

    Gezien een lijst met GTF's, die bij het samenvoegen opnieuw werden geschat, kunnen gebruikers het onderstaande protocol volgen om DESeq2 te gebruiken voor differentiële expressie-analyse. Om GTF's van onbewerkte reads te genereren, volgt u het StringTie-protocolpapier (tot aan de Ballgown-stap).


    Oefening 2: Differentiële expressie meten

    Selecteer de referentiesequentie opnieuw en u zou nu annotatietracks op expressieniveau moeten zien voor elke sampleconditie die erop is geladen. Als je deze tracks niet kunt zien, controleer dan of ze zijn ingeschakeld op het tabblad Annotatie en Tracks rechts van de sequentieviewer. Om genen te vinden die differentieel tot expressie worden gebracht tussen de twee monstercondities, ga naar: Annoteren en voorspellen → Expressieniveaus vergelijken.

    Controleer of de tracks die u in de vorige oefening hebt gemaakt ter vergelijking zijn gekozen: ze zouden moeten zijn Track 1 – Expressie: Sample_condition_1 en Track 2 – Expressie: Sample_condition_2. Kies om te vergelijken transcripties genormaliseerd door Mediaan van genexpressieverhoudingen – dit is de aanbevolen methode (zie metingen op expressieniveau voor meer informatie over deze methoden). Klik Oke.

    U ziet nu een derde annotatietrack toegevoegd aan uw referentiereeks, genaamd “Diff Expression: Sample_condition_1 vs Sample_condition_2”, en net als bij de individuele tracks op expressieniveau, worden annotaties gekleurd volgens de resultaten. Beweeg de muisaanwijzer over een van de annotaties en u ziet een pop-upvenster met de onbewerkte lees- en transcripttellingen voor elke monstervoorwaarde, plus een lijst met differentiële expressiescores.

    De Differentiële expressie p-waarde vertelt u of de waargenomen differentiële expressie statistisch significant is. De Differentiële expressie vertrouwen is de negatieve base 10 log van de p-waarde, aangepast om negatief te zijn voor genen die onder tot expressie worden gebracht in monster 2 in vergelijking met monster 1, of positief voor genen die tot overexpressie worden gebracht. Standaard gebruikt Geneious deze waarde om de annotaties in te kleuren: van blauw voor genen die onder tot expressie worden gebracht, tot wit voor genen die niet differentieel tot expressie worden gebracht, tot rood voor genen die tot overexpressie worden gebracht.

    Probeer te filteren op basis van dit veld om alle genen te vinden die een significant hogere expressie vertonen in monsterconditie 2 dan in monsterconditie 1: Geef de Annotaties en tracks tabblad rechts van de reeksviewer en typ “Differentiële expressie vertrouwen'8221>2 (inclusief de aanhalingstekens) in het zoekvak, zoals in de onderstaande schermafbeelding. Een betrouwbaarheidsniveau voor differentiële expressie boven 2 komt overeen met een p-waarde van minder dan 0,01. U zou nu slechts 29 annotaties moeten zien, die allemaal roze of rood van kleur zijn, wat wijst op een significant hogere expressie van dat gen in monsterconditie 2 dan in monsterconditie 1 (bij p<0.01). Om door deze annotaties te bladeren, klikt u op de pijlen rechts van de naam van de annotatietrack op het tabblad Annotaties en sporen.

    U kunt de resultaten ook in tabelvorm weergeven door op de annotaties tab boven de sequentieviewer en toont alleen de annotatietrack “Diff Expression'8221. Mogelijk moet u de kolommen voor Differential Expression Confidence en Ratio toevoegen door op de Kolommen knop. Om de resultaten te sorteren in de volgorde van meest overexpressie in voorbeeld 2 → meest onderexpressie in voorbeeld 2, klikt u tweemaal op de kop '&8220Differential Expression Confidence”, zodat de hoogste waarden bovenaan staan.

    U kunt deze tabel desgewenst in .csv-indeling exporteren door te klikken op Tabel exporteren (dit kan onder de dubbele pijlen >> staan).


    Opmerking: er zijn twee voorgestelde manieren om schattingen te importeren voor gebruik met differentiële genexpressie (DGE)-methoden. De eerste methode, die we hieronder laten zien voor edgeR en voor DESeq2, is om de geschatte tellingen op genniveau van de kwantificatietools te gebruiken, en bovendien om de overvloedschattingen op transcriptniveau te gebruiken om een ​​offset op genniveau te berekenen die corrigeert voor wijzigingen in de gemiddelde transcriptlengte over monsters. De onderstaande codevoorbeelden voeren deze stappen voor u uit, houden de juiste matrices bij en berekenen deze offsets. Voor edgeR je moet een matrix toewijzen aan y$offset , maar de functie DESeqDataSetFromTximport verzorgt de creatie van de offset voor u. Laten we deze methode "originele tellingen en offset”.

    De tweede methode is om het tximport-argument countsFromAbundance="lengthScaledTPM" of "scaledTPM" te gebruiken en vervolgens de telmatrix op genniveau txi$counts direct te gebruiken zoals u zou doen met een normale telmatrix met deze software. Laten we deze methode "bias gecorrigeerde tellingen zonder een offset

    Opmerking: Geef de oorspronkelijke tellingen op genniveau niet handmatig door aan downstream-methoden zonder een offset. Het enige geval waarin dit zinvol zou zijn, is als er geen lengteafwijking is voor de tellingen, zoals gebeurt in 3'-gelabelde RNA-seq-gegevens (zie de sectie hieronder). De oorspronkelijke tellingen op genniveau zijn in txi$counts wanneer tximport werd uitgevoerd met countsFromAbundance="no" . Dit is simpelweg het doorgeven van de gesommeerde geschatte transcripttellingen en corrigeert niet voor mogelijk differentieel isovormgebruik (de offset), wat het punt is van de tximport methoden (Soneson, Love en Robinson 2015) voor analyse op genniveau. Het doorgeven van ongecorrigeerde tellingen op genniveau zonder een offset wordt niet aanbevolen door de tximport pakket auteurs. De twee methoden die we hier aanbieden zijn: “originele tellingen en offset" of "bias gecorrigeerde tellingen zonder een offset”. Het doorgeven van txi aan DESeqDataSetFromTximport zoals hieronder beschreven is correct: de functie creëert de juiste offset voor u om differentiële expressie op genniveau uit te voeren.


    3'-gelabeld RNA-seq

    Als u 3'-gelabelde RNA-seq-gegevens hebt, zal het corrigeren van de tellingen voor genlengte een vertekening in uw analyse veroorzaken, omdat de tellingen geen lengtevertekening hebben. In plaats van de standaard pijplijn met volledige transcriptielengte te gebruiken, raden we aan om de originele tellingen te gebruiken, b.v. txi$counts als counts matrix, b.v. verstrekken aan DESeqDataSetFromMatrix of naar de edgeR of limma functies zonder een offset te berekenen en zonder gebruik te maken van teltVanOvervloed.


    Tellen normalisatie van Mov10 dataset met DESeq2

    Nu we de theorie van telling-normalisatie kennen, zullen we de tellingen voor de Mov10-dataset normaliseren met DESeq2. Dit vereist een paar stappen:

    1. Zorg ervoor dat de rijnamen van het metadata-dataframe aanwezig zijn en in dezelfde volgorde staan ​​als de kolomnamen van het counts-dataframe.
    2. Een DESeqDataSet-object maken
    3. Genereer de genormaliseerde tellingen

    1. Match de metadata en telt data

    We moeten er altijd voor zorgen dat we voorbeeldnamen hebben die overeenkomen tussen de twee bestanden, en dat de voorbeelden in dezelfde volgorde staan. DESeq2 geeft een foutmelding als dit niet het geval is.

    Als uw gegevens niet overeenkomen, kunt u de functie match() gebruiken om ze opnieuw te rangschikken.

    Stel dat we voorbeeldnamen hebben die overeenkomen in de tellingenmatrix en het metadatabestand, maar ze staan ​​in een andere volgorde. Schrijf de regel(s) code om een ​​nieuwe matrix te maken met kolommen die zo zijn geordend dat ze identiek zijn aan de rijnamen van de metadata.

    2. Creëer DESEq2-object

    Bioconductor-softwarepakketten definiëren en gebruiken vaak een aangepaste klasse binnen R voor het opslaan van gegevens (invoergegevens, tussentijdse gegevens en ook resultaten). Deze aangepaste gegevensstructuren zijn vergelijkbaar met lijsten omdat ze meerdere verschillende gegevenstypen/-structuren kunnen bevatten. Maar in tegenstelling tot lijsten, hebben ze vooraf gespecificeerd: dataslots, die specifieke soorten/klassen gegevens bevatten. De gegevens die in deze vooraf gespecificeerde slots zijn opgeslagen, kunnen worden geopend met behulp van specifieke pakketgedefinieerde functies.

    Laten we beginnen met het maken van het DESeqDataSet-object, en dan kunnen we wat meer praten over wat erin is opgeslagen. Om het object te maken, hebben we de . nodig telmatrix en de metagegevens tabel als invoer. We moeten ook een specificeren ontwerp formule:. De ontwerpformule specificeert de kolom(men) in de metadatatabel en hoe deze in de analyse moeten worden gebruikt. Voor onze dataset hebben we maar één kolom waarin we geïnteresseerd zijn, namelijk:

    monstertype . Deze kolom heeft drie factorniveaus, wat DESeq2 vertelt dat we voor elk gen genexpressieverandering willen evalueren met betrekking tot deze verschillende niveaus.

    Onze telmatrixinvoer wordt opgeslagen in het txi-lijstobject. We moeten dat dus specificeren met de functie DESeqDataSetFromTximport(), die de component counts extraheert en de waarden rondt op het dichtstbijzijnde gehele getal.

    OPMERKING: Omdat we in de vorige les een datavariabele hadden gemaakt die de tellingen bevat, hadden we die ook als invoer kunnen gebruiken. In dat geval zouden we echter de functie DESeqDataSetFromMatrix() willen gebruiken.

    U kunt DESeq-specifieke functies gebruiken om toegang te krijgen tot de verschillende slots en informatie op te halen. Stel dat we bijvoorbeeld de oorspronkelijke telmatrix willen ophalen, dan gebruiken we de functie counts() (Opmerking: we nesten het in de View()-functie zodat we het resultaat in de scripteditor kunnen bekijken in plaats van in de console) :

    Terwijl we de workflow doorlopen, zullen we relevante functies gebruiken om te controleren welke informatie in ons object is opgeslagen.

    3. Genereer de genormaliseerde Mov10-tellingen

    De volgende stap is om de telgegevens te normaliseren om eerlijke genvergelijkingen tussen monsters te maken.

    Om de . uit te voeren mediaan van ratio's methode van normalisatie, DESeq2 heeft een enkele schattingSizeFactors() functie die groottefactoren zal genereren. We zullen deze functie in het onderstaande voorbeeld demonstreren, maar in een typische RNA-seq-analyse wordt deze stap automatisch uitgevoerd door de functie DESeq(), die we later zullen bespreken.

    Door de resultaten weer toe te wijzen aan het dds-object, vullen we de slots van het DESeqDataSet-object met de juiste informatie. We kunnen de normalisatiefactoren van elk monster bekijken met behulp van:

    Om nu de genormaliseerde counts-matrix uit dds op te halen, gebruiken we de functie counts() en voegen we het argument normalized=TRUE toe.

    We kunnen deze genormaliseerde datamatrix opslaan in een bestand voor later gebruik:

    OPMERKING: DESeq2 gebruikt niet echt genormaliseerde tellingen, maar gebruikt de onbewerkte tellingen en modelleert de normalisatie binnen het gegeneraliseerde lineaire model (GLM). Deze genormaliseerde tellingen zullen nuttig zijn voor downstream-visualisatie van resultaten, maar kunnen niet worden gebruikt als invoer voor DESeq2 of andere tools die differentiële expressieanalyse uitvoeren die het negatieve binomiale model gebruiken.

    Deze les is ontwikkeld door leden van het docententeam van de Harvard Chan Bioinformatics Core (HBC). Dit zijn open access-materialen die worden gedistribueerd onder de voorwaarden van de Creative Commons Attribution-licentie (CC BY 4.0), die onbeperkt gebruik, distributie en reproductie in elk medium toestaat, op voorwaarde dat de oorspronkelijke auteur en bron worden vermeld.


    Bekijk de video: Contribution Margin explained (Januari- 2022).