Search for reads aligned to multiple chromosomes
1
0
Entering edit mode
5.0 years ago
oleglovky • 0

Hi,

I have sam/bam files containing rna-seq reads aligned to a genome. I need to find reads that are aligned to more than one chromosome. Is there a way (using samtools or other software) to achieve this?

EDIT: Even better would be a way to graphically view (similarly to IGV) aligned reads to more than one chromosome at a time and compare between them.

Thanks.

Assembly rna-seq genome alignment gene • 2.0k views
ADD COMMENT
0
Entering edit mode

First thoughts on that - I think you would be best off writing a python script to do it. Python has modules that allow parsing BAM/SAM files easy and then I expect you can get co ordinates out easily too. Then you can work out which reads map to multiple chromosomes. It would make a fun programming task.

ADD REPLY
0
Entering edit mode

Thanks for the tip! I have little experience with python, but I will surely look into this.

ADD REPLY
2
Entering edit mode
5.0 years ago
Macspider ★ 3.4k

I need to find reads that are aligned to more than one chromosome

If with this you mean chimeric reads, then the easiest yet time consuming step would be to re-run the mapping (perhaps with TopHat-fusion). Fusion related options:

--fusion-search
--fusion-anchor-length         <int>       [ default: 20               ]
--fusion-min-dist              <int>       [ default: 10000000         ]
--fusion-read-mismatches       <int>       [ default: 2                ]
--fusion-multireads            <int>       [ default: 2                ]
--fusion-multipairs            <int>       [ default: 2                ]
--fusion-ignore-chromosomes    <list>      [ e.g, <chrM,chrX>          ]

If, instead, you mean pairs that have mates mapping in different scaffolds, then it's just a matter of setting the right bitwise flag on samtools view. That would be:

samtools view -F 0x2 -F 0x0100 -F0x0400 -F0x4 -F 0x8

(There must be a simpler way to write this flag selector but I came to this one recently by myself and seems to work at least for my purposes)

ADD COMMENT

Login before adding your answer.

Traffic: 3251 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6