Vraag:
Hoe cross alignments uit een BED-bestand filteren?
SmallChess
2017-05-19 10:49:47 UTC
view on stackexchange narkive permalink

Ik heb een BAM-bestand:

  @SQ SN: chr1 LN: 248956422 @ SQ SN: chrx LN: 248956423ST-E00110: 348: HGVKKALXX: 1: 1201: 5822: 48670 323 Chr1 9999 0 67H66M16H chrx 1000 0 GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ NM: i: 0 MD: Z: 66 AS: i: 66 XS: i: 65 SA: Z: chr5,18606834, -, 73S76M, 34,0; RG: Z: g1  

Er is een read uitgelijnd op chr1 , en het is mat uitgelijnd op chrx .

Ik heb een BED-bestand:

  chr1 0 100000 TestOnly  

Ik wil graag alles eruit filteren wat buiten mijn BED valt regio, die kruisuitlijningen bevat. In mijn voorbeeld, hoewel mijn read is uitgelijnd op chr1 maar het is mate niet. Ik wil dit niet lezen.

Wanneer ik dat doe:

samtools view -L test.bed test.bam

het commando geeft me de read omdat het de cross-alignments niet controleert.

Mijn oplossing:

samtools view -L test.bed test.bam | grep -v chrx

maar dit is erg traag en onhandig. In mijn productiepijplijn zou ik zoiets moeten doen als:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

V: Is er een betere oplossing?

Een antwoord:
#1
+6
terdon
2017-05-19 22:44:29 UTC
view on stackexchange narkive permalink

Volgens de SAM-specificatie is het derde veld van een SAM-regel ( RNAME ):

RNAME: Referentiereeks NAME van de uitlijning. Als @SQ-headerregels aanwezig zijn, moet RNAME (indien niet ‘*’) aanwezig zijn in een van de SQ-SN-tags. Een niet-toegewezen segment zonder coördinaat heeft een ‘*’ bij dit veld. Een niet-toegewezen segment kan echter ook een gewone coördinaat hebben, zodat het na het sorteren op een gewenste positie kan worden geplaatst. Als RNAME '*' is, kunnen er geen aannames worden gedaan over POS en CIGAR.

En het zevende veld is (nadruk van mij, ontbreekt "aan" die van hen):

RNEXT: Referentiereeksnaam van de primaire uitlijning van de NEXT die in het sjabloon wordt gelezen. Voor het laatst gelezen, is het volgende dat het eerst wordt gelezen in de sjabloon. Als @SQ-headerregels aanwezig zijn, moet RNEXT (indien niet ‘*’ of ‘=’) aanwezig zijn in een van de SQ-SN-tags. Dit veld wordt ingesteld als ‘*’ als de informatie niet beschikbaar is, en wordt ingesteld als ‘=’ als RNEXT identiek is aan RNAME . Als niet ‘=’ en de volgende gelezen in de sjabloon één primaire mapping heeft (zie ook bit 0x100 in FLAG), is dit veld identiek aan RNAME op de primaire regel van de volgende gelezen. Als RNEXT '*' is, kunnen er geen aannames worden gedaan over PNEXT en bit 0x20

Dus je wilt die regels verwijderen waarvan het 7e veld niet = is en, voor het geval dat, die regels waarvan het zevende veld niet = en is, zijn niet hetzelfde als het derde veld. U kunt daarom zoiets als dit gebruiken:

  samtools view -L test.bed test.bam | awk '$ 7 == "=" || $ 3 == $ 7  

En, om weer op te slaan als een bam-bestand:

  samtools bekijk -L test.bed test.bam | awk '$ 7 == "=" && $ 3 == $ 7 | samtolls-weergave -b > fixed.bam  

Op een aparte opmerking: het is zeer zelden nodig om meerdere grep-commando's op deze manier aan elkaar te koppelen. U kunt \ | (of | met de -E of -P -opties) gebruiken om ze te scheiden. Iets als:

  samtools view -L test.bed test.bam | grep -v 'chrx \ | chr2 \ | chr10 \ | chrN'  

Of

  samtools-weergave -L test.bed test.bam | grep -Ev 'chrx | chr2 | chr10 | chrN'  
Als je het op deze manier doet, mist het bestand `fixed.bam` de header, wat naar mijn ervaring veel problemen oplevert. Ik raad aan om de koptekst altijd weer toe te voegen; ofwel door `-h` op te geven bij het lezen van de originele BAM, of door het afzonderlijk toe te voegen:` (samtools view -H infile.bam; samtools view…)> samtools view -b> outfile.bam`.


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