Otázka:
Jak mohu zrychlit volání / opravy INDEL na souborech BAM?
gringer
2017-06-05 04:25:12 UTC
view on stackexchange narkive permalink

Příkaz samtools mpileup má celkem pěknou vlastnost, že je schopen opravit chyby mapování spojené se špatným zarovnáním INDELů. Ve výchozím nastavení nebude příkaz mpileup fungovat pro čtení, která mají více než 250násobné pokrytí referenčního genomu. I když lze tento limit zvýšit, velmi vysoké pokrytí způsobí, že se program mpileup zastaví, takže by bylo hezké vědět, jestli existuje nějaký snadný způsob, jak to zrychlit.

Chcete-li přidat trochu více V tomto kontextu jsem to dělal s čteními mitochondriálního genomu, které byly extrahovány jak z celogenomového sekvenování Illumina (pokrytí ~ 1000X), tak z cíleného amplikonového sekvenování provedeného na IonTorrent (pokrytí až ~ 4000X).

Vidím, že @rightskewed zmínil schopnost převzorkování samtoolů s samtools view -s <float> (viz zde), což vypadá, že by to mohlo fungovat jako řešení pokud je použit před operací mpileup.

Jeden odpovědět:
gringer
2017-06-05 04:34:22 UTC
view on stackexchange narkive permalink

Když jsem měl tento problém před několika lety, nevěděl jsem o podvzorkování samtoolů, takže jsem nakonec napsal vlastní metodu digitální normalizace, abych se vypořádal s mapovanými čteními. Tato metoda snižuje pokrytí genomu, ale zachovává čtení tam, kde je nízké pokrytí.

Protože jsem pracoval s čteními IonTorrent (která mají proměnnou délku), přišel jsem s nápadem vybrat nejdelší čtení, které je mapováno na každé umístění v genomu (za předpokladu, že takové čtení existovalo). To znamenalo, že vysoce variabilní pokrytí pro různé vzorky (někdy až 200x, někdy až 4000X) bylo vyrovnáno na mnohem konzistentnější pokrytí kolem 100-200X. Tady je jádro kódu Perlu, který jsem napsal:

  if (($ F [2] ne $ seqName) || ($ F [3]! = $ Pos) || (délka ($ bestSeq) < = length ($ F [9]))) {if (length ($ bestSeq) == length ($ F [9])) {## vzorkování nádrže s velikostí nádrže 1 ## viz https : //en.wikipedia.org/wiki/Reservoir_sampling ## * s pravděpodobností 1 / i, ponechat novou položku namísto aktuální položky $ seenCount ++; if (! rand ($ seenCount)) {## tj. pokud rand ($ seenCount) == 0, pak pokračujte s výměnou další; }} else {$ seenCount = 1; } if (($ F [2] ne $ seqName) || ($ F [3]! = $ pos)) {if ($ output eq "fastq") {printSeq ($ bestID, $ bestSeq, $ bestQual); } elsif ($ output eq "sam") {print ($ bestLine); }} $ seqName = $ F [2]; $ pos = $ F [3]; $ bestLine = $ line; $ bestID = $ F [0]; $ bestFlags = $ F [1]; $ bestSeq = $ F [9]; $ bestQual = $ F [10];  

Celý kód (který funguje jako filtr pro nekomprimované soubory SAM) najdete zde.

K objasnění: toto je kód Perlu?
ano, Perl kód. Můžu to poznat podle symbolů dolaru všude, ale myslím, že to není zřejmé pro lidi, kteří Perla neznají.
Mluvím sice Perl, ale je to stěží jediný jazyk s dolary před identifikátory. ;-)
Vskutku. Moje mysl byla prázdná, ale moje žena mi připomínala bash / sh.


Tyto otázky a odpovědi byly automaticky přeloženy z anglického jazyka.Původní obsah je k dispozici na webu stackexchange, za který děkujeme za licenci cc by-sa 3.0, pod kterou je distribuován.
Loading...