Vraag:
Is er een gemakkelijke manier om een ​​samenvatting te maken van een VCF-bestand (v4.1) met structurele variaties?
Kamil S Jaron
2017-05-29 02:40:14 UTC
view on stackexchange narkive permalink

Ik heb een aantal vcf-bestanden (v4.1) met structurele variaties van een aantal niet-modelorganismen (d.w.z. er zijn geen bekende varianten). Ik ontdekte dat er nogal wat tools zijn om vcf-bestanden te manipuleren, zoals VCFtools, R-pakket vcfR of python-bibliotheek PyVCF. Geen van hen lijkt echter een korte samenvatting te geven, zoiets als (liefst ook gecategoriseerd op grootte):

  type countDEL xINS yINV z ....  

Is er een tool of een functie die ik over het hoofd heb gezien die samenvattingen van deze stijl oplevert?

Ik weet dat het vcf-bestand slechts een gewoon tekstbestand is en of ik REF en ALT -kolommen Ik zou een script moeten kunnen schrijven dat de klus zal klaren, maar ik hoopte dat ik kon vermijden om mijn eigen parser te schrijven.

--- bewerken ---

Tot dusver lijkt het erop dat alleen tool die tot doel heeft samenvattingen te maken (@gringer antwoord) niet werkt op vcf v4.1. Andere tools zouden slechts een gedeeltelijke oplossing bieden door een bepaald varianttype te filteren. Daarom accepteer ik mijn eigen parser perl / R-oplossingen, totdat er een werkende tool zal zijn voor statistieken van vcf met structurele varianten .

Wat voor soort samenvatting zoekt u? Alleen ruwe tellingen (van het aantal invoegingen / verwijderingen)? Ik ben gewoon nieuwsgierig, omdat ik persoonlijk denk dat ik geïnteresseerd zou zijn om te weten hoe die variatie zich ruimtelijk langs een genoom / contig / sequentie verspreidt, zou dat interessant voor je zijn?
Ik vergelijk een paar relatieve soorten. Ik heb echter nog geen homologe regio's geïdentificeerd (de annotatie is tot dusver in aanbouw), daarom is de positionering van SV's op dit moment niet erg relevant. Wat ik echter al kan doen, is kijken of er op zijn minst een verschil is in tellingen / grootte van SV's tussen soorten.
Misschien zelfs beter dan samenvatting zou zijn om het type SV en de grootte voor elke SV-oproep te krijgen. Dan zou ik een samenvatting in R triviaal kunnen zijn.
Je hebt dit waarschijnlijk al gezien, maar `vcftools` heeft een broer of zus genaamd` bcftools`, die [een zoekfunctie heeft, waarmee gebruikers een VCF / BCF kunnen opvragen om velden en informatie eruit te halen en hun eigen formaat uit te voeren] (https : //samtools.github.io/bcftools/bcftools-man.html#query). Het doet misschien niet precies wat je wilt, maar het kan je wel heel dichtbij brengen (genoeg om misschien wat nabewerking in R nodig te hebben?).
Geweldig, dat klinkt veelbelovend. Ik zal er naar kijken.
Vijf antwoorden:
#1
+8
gringer
2017-05-29 10:16:36 UTC
view on stackexchange narkive permalink

Volgens de bcftools man-pagina, is het in staat om statistieken te produceren met het commando bcftools stats . Als u dit zelf uitvoert, zien de statistieken eruit zoals u vraagt:

  # Dit bestand is gemaakt door bcftools stats (1.2-187-g1a55e45 + htslib-1.2.1-256-ga356746) en kan worden geplot met plot-vcfstats. # De opdrachtregel was: bcftools stats OVLNormalised_STARout_KCCG_called.vcf.gz ## Definitie van sets: # ID [2] id [3] door tabs gescheiden bestandsnamenID 0 OVLNormalised_STARout_KCCG_called.vcf.gz # SN , Samenvattingsnummers: # SN [2] id [3] sleutel [4] waardeSN 0 aantal samples: 108SN 0 aantal records: 333SN 0 aantal no-ALT's: 0SN 0 aantal SNP's: 313SN 0 aantal MNP's: 0SN 0 aantal indels: 20SN 0 aantal andere: 0SN 0 aantal multi-parallelle sites: 0SN 0 aantal multi-parallelle SNP-sites: 0 # TSTV, overgangen / transversies: # TSTV [2] id [3] ts [4] tv [5 ] ts / tv [6] ts (1e ALT) [7] tv (1e ALT) [8] ts / tv (1e ALT) TSTV 0302 11 27.45 302 11 27.45 # SiS, Singleton stats: ... # IDD, InDel-distributie: # IDD [2] id [3] lengte (deleti ons negatief) [4] countIDD 0-9 1IDD 0 -2 4IDD 0-1 6IDD 0 1 4IDD 0 2 1IDD 0 3 3IDD 0 4 1 # ST, vervangingstypes: # ST [2] id [3] type [4] countST 0 A>C 2ST 0 A>G 78ST 0 A>T 2ST 0 C>A 5ST 0 C>G 0ST 0 C>T 66ST 0 G>A 67ST 0 DPGiliqoxvy # 0ST 0 G>C 0ST 0 G>C 0ST 0 G>A 67ST 0 G>C 0ST 0 G>A ] bin [4] aantal genotypen [5] fractie genotypen (%) [6] aantal sites [7] fractie van sites (%) DP 0 >500 0 0,000000 333100,000000  
Ik doet bijna precies wat ik zou willen zien. Maar het is blijkbaar niet gemaakt voor SNV's, SV's. Het heeft breekpunten verkeerd gelabeld als indels.
Het annoteert wat expliciet in het VCF-bestand staat, en dat zijn SNV's (en INDEL's). Als je een structurele variantanalyse wilt (d.w.z. op grotere schaal dan enkele nucleotiden), dan moet je iets gebruiken dat meer doet dan een samenvatting van het VCF-bestand. Inversies, grootschalige verwijderingen en breekpunten maken geen deel uit van het "standaard" VCF-formaat.
Ze maken deel uit van 4.1. Het wordt beschreven in de documentatie: https://samtools.github.io/hts-specs/VCFv4.1.pdf
Juist, ik begrijp het. Ik denk dat het dan een bcftools-bug is: https: //github.com/samtools/bcftools/issues
Ok, bevestigd: bcftools ondersteunt GEEN structurele varianten https://github.com/samtools/bcftools/issues/623
#2
+6
Kamil S Jaron
2017-05-29 11:22:59 UTC
view on stackexchange narkive permalink

De "mijn eigen parser" -oplossingen. De informatie waar ik naar zocht in een deel van de kolom INFO , namelijk in variabelen SVLEN en SVTYPE.

zeer snelle SV types + counts (door @ user172818 in commnent):

  zcat var.vcf.gz | perl -ne 'print "$ 1 \ n" if / [; \ t] SVTYPE = ([^; \ t] +) /' | sorteren | uniq -c  

vrij trage SV-typen + counts + formaten:

  SV_colnames <- c ('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1') ssplit <- functie (s, split = '=') {unlist (strsplit (s, split = split))} # opmerking, hoofdletters respecteren gewoon de originele naamgevingsconventies van de VCF filegetSVTYPE <- functie (info) {ssplit (grep ("SVTYPE", info, waarde = T)) [2]} getSVLEN <- functie (info ) {SVLEN <- ssplit (grep ("SVLEN", info, waarde = T)) ifelse (lengte (SVLEN) == 0, NA, as.numeric (SVLEN [2]))} load_sv <- functie (bestand) {vcf_sv_table <- read.table (file, stringsAsFactors = F) colnames (vcf_sv_table) <- SV_colnames # mogelijke filtering # vcf_sv_table <- vcf_sv_table [vcf_sv_table $ FILTER == 'PASS',] vcf_sv_table_info <- strsplit (vcf_sv_table $ INFO, ' ; ') vcf_sv_table $ SVTYPE <- unlist (lapply (vcf_sv_table_info, getSVTYP E)) vcf_sv_table $ SVLEN <- unlist (lapply (vcf_sv_table_info, getSVLEN)) return (vcf_sv_table)} sv_df <- load_sv ('my_sv_calls./TvcYf')
Als je SVTYPE alleen wilt tellen, zal deze opdrachtregel sneller zijn: `zcat var.vcf.gz | perl -ne 'print "$ 1 \ n" if / [; \ t] SVTYPE = ([^; \ t] +) /' | sorteren | uniq -c`. SV-bestanden zijn klein, dus je Rscript is prima, maar over het algemeen is R erg traag voor tekstverwerking. Perl / Python / etc heeft de voorkeur als je te maken hebt met enorme VCF's.
De R-oplossing laat deuren open voor extra spelen met lengte / positie / sv-type, maar je hebt absoluut gelijk over de uitvoering.
#3
+2
eastafri
2017-05-29 10:23:32 UTC
view on stackexchange narkive permalink

Misschien wilt u Bio-VCF proberen. Uit de beschrijving van de auteur

Bio-vcf is een nieuwe generatie VCF-parser, filter en converter. Bio-vcf is niet alleen erg snel voor genoombrede (WGS) data, het wordt ook geleverd met een erg goede filter-, evaluatie- en herschrijftaal en het kan elk type tekstuele data uitvoeren, inclusief VCF-header en inhoud in RDF en JSON.

Bovendien is het een zeer snelle parser. De herschrijf DSL kan u helpen uw filtering en behoeften aan te passen.

Heeft Bio-vcf enige functionaliteit voor structurele variatie zoals gevraagd door het OP? Een snelle scan van de GitHub-pagina onthult niets. Hoewel veel van deze tools nuttig zijn voor het filteren van SNV's en korte indels, doen de meeste niets van waarde voor SV's.
#4
+2
Medhat Helmy
2019-09-17 01:19:02 UTC
view on stackexchange narkive permalink

Voor structurele varianten denk ik dat je SURVIVOR zoals SURVIVOR-statistieken kunt gebruiken die specifiek voor een dergelijk doel zijn ontworpen (statistieken over SV-bestand).

#5
+1
Panwen Wang
2019-10-22 04:01:05 UTC
view on stackexchange narkive permalink

Stel dat je een INFO-veld hebt met de naam SVTYPE om het type structuurvariant aan te geven.

Hier is hoe je de statistieken eenvoudig kunt opvragen:

  1. Installeer vcfstats : pip install vcfstats
  2. Definieer een macro om de SVTYPE -informatie te extraheren, zeg in svtype.py :
  van vcfstats.macros import cat @ catdef SVTYPE (variant): return variant.INFO [ 'SVTYPE']  
  1. Genereer de statistieken:
  vcfstats --vcf <your vcf file> \ --outdir <the output directory> \ --macro svtype.py \ --formule 'COUNT (1) ~ SVTYPE [DEL, INS, INV ...]' \ # gewenste svtypes - titel 'Aantal SV-typen'  

Vervolgens heeft u in de uitvoermap het statistiekenbestand met de naam Number_of_SV_Types.txt en ook een plot ervoor.

Controleer de details op https: // git hub.com/pwwang/vcfstats

Het voordeel van deze oplossing is dat er een standaard parser wordt gebruikt. Het nadeel is dat het verschillende stappen vereist (in tegenstelling tot `zcat var.vcf.gz | perl -ne 'print" $ 1 \ n "if / [; \ t] SVTYPE = ([^; \ t] +) /' | sort | uniq -c` die de samenvatting in één keer zal maken).


Deze Q&A is automatisch vertaald vanuit de Engelse taal.De originele inhoud is beschikbaar op stackexchange, waarvoor we bedanken voor de cc by-sa 3.0-licentie waaronder het wordt gedistribueerd.
Loading...