Saving to BAM files from a BioMap object ?

1 view (last 30 days)
Hi,
I noticed that in the Release Notes of the 2013a bioinformatics toolbox it is mentioned that now it is possible to write BAM files, but I couldn't find any example. Could someone give me a simple example of using the new feature. For example if I have a BAM file with the ScannedDictionary 'chrI' and 'chrII', how can I split this into 2 smaller BAM files?
Thanks,
Razvan

Accepted Answer

Lucio Cetto
Lucio Cetto on 6 May 2013
You can not only split, but you can use this convenient method to filter (or preprocess) your files, for example you can filter by quality or by flags very easily (subsetting a file indexed BioMap without the need to create temporal files). In the next example I load a 2.6Gb BAM file and I save only reads which lengths are 75 and were mapped to chromosome 1.
>> bhe1 = BioMap('E:\rnaseq\humanembryonic\SRX026674\SRR065504\accepted_hits.bam')
bhe1 =
BioMap with properties:
SequenceDictionary: {1x25 cell}
Reference: [41927566x1 File indexed property]
Signature: [41927566x1 File indexed property]
Start: [41927566x1 File indexed property]
MappingQuality: [41927566x1 File indexed property]
Flag: [41927566x1 File indexed property]
MatePosition: [41927566x1 File indexed property]
Quality: [41927566x1 File indexed property]
Sequence: [41927566x1 File indexed property]
Header: [41927566x1 File indexed property]
NSeqs: 41927566
Name: ''
>> len = bhe1.getStop-bhe1.getStart+1; % using short reads mapped to transcripts
bnhe1 = getSubset(bhe1,len==75)
bnhe1 =
BioMap with properties:
SequenceDictionary: {1x25 cell}
Reference: [30558482x1 File indexed property]
Signature: [30558482x1 File indexed property]
Start: [30558482x1 File indexed property]
MappingQuality: [30558482x1 File indexed property]
Flag: [30558482x1 File indexed property]
MatePosition: [30558482x1 File indexed property]
Quality: [30558482x1 File indexed property]
Sequence: [30558482x1 File indexed property]
Header: [30558482x1 File indexed property]
NSeqs: 30558482
Name: ''
>> bnhe2 = getSubset(bnhe1,'SelectReference',1)
bnhe2 =
BioMap with properties:
SequenceDictionary: '1'
Reference: [2907478x1 File indexed property]
Signature: [2907478x1 File indexed property]
Start: [2907478x1 File indexed property]
MappingQuality: [2907478x1 File indexed property]
Flag: [2907478x1 File indexed property]
MatePosition: [2907478x1 File indexed property]
Quality: [2907478x1 File indexed property]
Sequence: [2907478x1 File indexed property]
Header: [2907478x1 File indexed property]
NSeqs: 2907478
Name: ''
>> write(bnhe2,'newbam')
>> ! ls -l E:\rnaseq\humanembryonic\SRX026674\SRR065504\new*
-rwxrwxrwa 183729053 Apr 9 20:38 E:/rnaseq/humanembryonic/SRX026674/SRR065504/newbam.bam
>> nb = BioMap('E:\rnaseq\humanembryonic\SRX026674\SRR065504\newbam.bam')
nb =
BioMap with properties:
SequenceDictionary: {1x25 cell}
Reference: [2907478x1 File indexed property]
Signature: [2907478x1 File indexed property]
Start: [2907478x1 File indexed property]
MappingQuality: [2907478x1 File indexed property]
Flag: [2907478x1 File indexed property]
MatePosition: [2907478x1 File indexed property]
Quality: [2907478x1 File indexed property]
Sequence: [2907478x1 File indexed property]
Header: [2907478x1 File indexed property]
NSeqs: 2907478
Name: ''
>> ! ls -l E:\rnaseq\humanembryonic\SRX026674\SRR065504\new*
-rwxrwxrwa 183729053 Apr 9 20:38 E:/rnaseq/humanembryonic/SRX026674/SRR065504/newbam.bam
-rwxrwxrwa 342176 Apr 9 20:39 E:/rnaseq/humanembryonic/SRX026674/SRR065504/newbam.bam.bai
-rwxrwxrwa 6621666 Apr 9 20:40 E:/rnaseq/humanembryonic/SRX026674/SRR065504/newbam.bam.linearindex
>>

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!