Vraag:
Een BAM-bestand van de ene referentie naar de andere converteren?
morgantaschuk
2017-05-17 01:00:43 UTC
view on stackexchange narkive permalink

Ik heb een set BAM-bestanden die zijn uitgelijnd met behulp van de NCBI GRCh37 menselijke genoomreferentie (met de chromosoomnamen als NC_000001.10), maar ik wil het analyseren met een BED-bestand dat de UCSC hg19-chromosoomnamen heeft (bijv. chr1 ). Ik wil bedtools gebruiken om alle on-target en off-target reads op te halen.

  1. Zijn NCBI en UCSC direct vergelijkbaar? Of moet ik de BAM / lift-over het BED opnieuw uitlijnen met de UCSC-referentie?
  2. Moet ik het BED-bestand of het BAM-bestand converteren? Iedereen hier gebruikt de UCSC-chromosoomnamen / -posities, dus ik moet de uiteindelijke bestanden toch naar UCSC converteren.
Drie antwoorden:
#1
+23
Devon Ryan
2017-05-17 02:04:43 UTC
view on stackexchange narkive permalink

Je bent de tweede persoon die ik ooit heb gezien met NCBI-"chromosoomnamen" (ze lijken meer op supercontig-ID's). Normaal gesproken zou ik u wijzen op een bron die toewijzingen tussen chromosoomnamen levert, maar aangezien niemand NCBI-namen heeft toegevoegd (maar misschien zal ik ze nu toevoegen), heeft u daar momenteel pech.

Hoe dan ook, de snelste manier om te doen wat je wilt, is door samtools view -H foo.bam > header te bekijken om de BAM-header te krijgen en vervolgens elke NCBI-"chromosoomnaam" te wijzigen in zijn overeenkomstige UCSC-chromosoomnaam. HERORDER DE LIJNEN NIET! Je kunt dan samtools reheader gebruiken en je bent klaar.

Waarom zou je je kunnen afvragen, zou dit werken? Het antwoord is dat chromosoom- / contig-namen in BAM-bestanden niet in elke uitlijning worden opgeslagen. In plaats daarvan worden de namen opgeslagen in een lijst in de koptekst en elke uitlijning bevat alleen de integer-index in die lijst (leesgroep-ID's zijn vergelijkbaar, voor wat het waard is). Dit leidt ook tot de waarschuwing hierboven tegen het opnieuw rangschikken van items, aangezien dat een ZEER gemakkelijke manier is om uitlijningen tussen chromosomen uit te wisselen.

Even terzijde, je zou er goed aan doen om over te schakelen naar gencode- of ensembl-chromosoomnamen, ze zijn iets meer samenhangend dan de iets_willekeurige puinhoop die aanwezig is in hg19 van UCSC.

Update : Omdat ik aardig ben, hier is de conversie tussen NCBI en UCSC. Merk op dat als u uitlijningen heeft met patches, er eenvoudig geen UCSC-equivalent is. Een van de vele redenen om UCSC niet te gebruiken (vermijd ook hun annotaties).

Hoe goed werkt dit? Heeft u benchmarks uitgevoerd? Vraag ik omdat ik heb geprobeerd om verschillende bedbestanden om te zetten tussen hg- en GRC-genomen en de drie tools die ik gebruikte, gaven allemaal heel verschillende resultaten. Dit soort mapping zou eigenlijk heel eenvoudig moeten zijn, maar het lijkt helemaal niet zo eenvoudig te zijn.
Voor gevallen waarin het slechts een naamswijziging is (de meeste gevallen), valt er niets te benchmarken. Voor gevallen waarin u bovendien positiewijzigingen heeft, heeft u een andere bron nodig (namelijk liftOver of crossmap).
Ja, het zijn tolgelden zoals liftOver en crossmap die ik heb gebruikt en waarmee ik problemen heb gevonden. Ik verwachtte dat dit een opgelost probleem zou zijn, maar elk van de drie tools die ik gebruikte, gaf helaas verschillende resultaten. Dat maakt me wantrouwend om de resultaten te gebruiken.
De resultaten zijn deterministisch, ze zouden hetzelfde moeten zijn, ongeacht de tool, op voorwaarde dat je dezelfde instellingen gebruikt.
Je zou het denken, ja. Maar dat zijn ze niet. Ik zou hier graag over willen chatten ([chatten]) (dit is een van de dingen die ik op het werk moet regelen) als je wilt, en ik zou heel blij zijn te horen dat ik het gewoon verkeerd deed, maar mijn voorlopige tests en verschillende discussies die ik online heb gelezen (zie [hier] (https://www.biostars.org/p/14187/#91490) bijvoorbeeld) suggereren dat het niet zo eenvoudig is als je zou denken.
Ik vermoed dat dit te wijten zal zijn aan tekortkomingen in de kettingbestanden, maar dat hoor ik graag.
#2
+3
chrisamiller
2017-05-17 23:56:14 UTC
view on stackexchange narkive permalink

De "juiste" oplossing zou herschikking zijn, maar dat is duur en de meesten van ons zouden die weg niet gaan. Mijn voorkeursoplossing zou zijn om het bedbestand te converteren, in tegenstelling tot de bam. Dit is waarom:

1) Herheadering van de bam betekent dat je reads hebt uitgelijnd met contigs zonder een corresponderende invoer in UCSC (zie Devon's lijst voor de mappings). Dit is een probleem omdat:

  • Sommige van die reads zouden waarschijnlijk ergens anders in kaart zijn gebracht als een referentie zonder deze contigs werd gebruikt.
  • Ik weet niet eens zeker wat er met die reads gebeurt nadat ze opnieuw zijn verzonden - ik denk dat ze moeten worden gemarkeerd als niet-toegewezen? Veel potentie voor smerigheid daar.

2) Het lijkt schoner om het bedbestand om te zetten van UCSC-> NCBI, waar je gegarandeerd bent dat elke ingang een "thuis" heeft. Daarna, nadat u uw informatie uit de bam heeft gehaald, kunt u de chromosoomnamen altijd terug converteren als dat nodig is.

Uitlijningen die uiteindelijk hun chromosoom verliezen, worden niet in kaart gebracht, hoewel ze niet de SAM niet-toegewezen vlag krijgen (in plaats daarvan worden ze weergegeven als zijnde op het '*' chromosoom). Maar ja, `samtools reheader` moet met uiterste voorzichtigheid worden gebruikt.
Verschillende tools zullen wiebelig zijn als u leest op `*` doorgeeft die niet als niet-toegewezen zijn gemarkeerd. `pysam` is er één. Ik weet dit uit bittere (bugfixing) ervaring.
#3
  0
Ian Sudbery
2017-08-08 16:02:03 UTC
view on stackexchange narkive permalink

Om erop te wijzen dat als je het antwoord van @Devon Ryan voor een ander organisme / assemblage wilt volgen, dat niet in zijn zeer bruikbare gekoppelde bron staat, je NCBI kunt downloaden naar UCSC contig to chromosome number mappings van https : //www.ncbi.nlm.nih.gov/assembly.

Ga naar de site en zoek naar je assembly-naam. Onderaan de pagina staat een vak genaamd "Globale assembly definition" met daarin een link getiteld "Download full sequence report".

Het gedownloade bestand bevat een tabel met:

  • Chromosoomnummers in "Sequence-Name" / "Assigned-molecule"
  • NCBI-namen in Refseq-Accn
  • UCSC-contignaam in "UCSC-style-name"


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...