Re: [bedtools-discuss] Not covered regions between bam and bed.

[Aaron Quinlan]

Hi Guillermo,

I need a bit kore information to make sense of this.  Could you please a) provide the exact commands you used in this analysis and b) provide a printout of the SAM records for the read ids in your example (e.g. samtools view aln.bam | grep "M01167:3:000000000-A37T8:1:1110:11018:11790/1")?

Thanks,

- Aaron

quinlanlab.org

On May 6, 2013, at 11:41 AM, gmarco <guim...@gmail.com> wrote:

Hello,

My input files are: 

• input.bam (paired end) Illumina alignment obtained with BWA.
  
• regions.bed file which contains the intervals of the current design

Here's a sample of regions.bed

chr1 10292367 10292512

chr1 10316285 10316401

chr1 10318531 10318750

chr1 10321943 10322048

chr1 10327418 10327636

chr1 10328190 10328341

The objective is to obtain which intervals of "regions.bed" design are not being covered in bam.

For that purpose I've tried converting the bam input into a bed file. Then intersecting both bed files with -v option.

The problem comes when I get this weird output in resulting bed file:

That 0 value is really weird...

chr1    10833   10931   M01167:3:000000000-A37T8:1:1110:11018:11790/1   0       +       10833   10931   0,0,0   1       98,     0,

chr1    10977   11027   M01167:3:000000000-A37T8:1:1101:10241:4209/2    0       +       10977   11027   0,0,0   1       50,     0,

And lines being repeated even if it's not the same read:

chr1    10292578        10292679        M01167:3:000000000-A37T8:1:2104:17089:16810/1   60      -       10292578        10292679        0,0,0   1       101,    0,

chr1    10292579        10292680        M01167:3:000000000-A37T8:1:2105:8403:22087/1    60      -       10292579        10292680        0,0,0   1       101,    0,

chr1    10292579        10292680        M01167:3:000000000-A37T8:1:2112:16677:2861/1    60      -       10292579        10292680        0,0,0   1       101,    0,

Thank you.

Best regards,

Guillermo.

--
You received this message because you are subscribed to the Google Groups "bedtools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bedtools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Aaron Quinlan]

Hi Guillermo,

As I understand your email, what you are looking for is the "sub"-intervals found in your target_regions.bed file that have no coverage whatsoever (e.g., 1-2, 8-10 in your example).  Is that correct?

- Aaron

On May 7, 2013, at 6:47 AM, Guillermo Marco <guim...@gmail.com> wrote:

Hi,

After realizing that I cannot really ensure if what I obtain is correct or not I made a small test setup:

1) short_alignment.bed (few lines converted from original BAM align):

[CODE]

chr1 10833 10931 M01167:3:000000000-A37T8:1:1110:11018:11790/1 0 +

chr1 10833 10883 M01167:3:000000000-A37T8:1:1101:10241:4209/2 0 +

chr1 10833 10933 M01167:3:000000000-A37T8:1:1113:11938:2404/2 16 +

chr1 10836 10937 M01167:3:000000000-A37T8:1:2109:28144:20710/1 37 +

chr1 11062 11154 M01167:3:000000000-A37T8:1:1101:10241:4209/1 0 -

[/CODE]

2) I manually edited my target_regions.bed for test purpose:

chr1    10130   11154

3) Launched the following commands:

intersectBed -a target_regions.bed -b short_alignment.bed -v

No output.

intersectBed -a target_regions.bed -b short_alignment.bed -v -f 1

chr1 10130 11154

With -f 1 option I can obtain intervals not totally being covered by 1 bp gap.

Is great however I want to obtain the part of the intervals not being covered.

Conclusion if my target region is from 1 to 10 and I map from 3 to 7 and would like to obtain intervals 1-2, 8-10.

I don't know if i'm explaining myself clearly.