Ideas on how to count the number of specific amplicons in aligned amplicon data

[swNGS]

I have a Haloplex library prep'd data. I would like to know if anyone can think of a way of counting the occurrence of a particular fragment that is derived from the original library once I've aligned the paired-end reads to the reference genome. I've assumed I can only do this from the aligned bam since the insert size, and hence the original fragment can be inferred his way.

So, an example as follows:
Fragment A drives from genomic position chr12 12345678 - 12345780

How many times does this feature in the aligned bam file?
I dont want to know the read depth across this fragment/insert ( could easily do that with bedtools).

The flip side of this question, is how about generating a frequency of occurrence of all inferred inserts in the aligned file? With genomic aligned coordrs of each inferred insert? Not a lookup this time, just a count of how many exist in the aligned bam file

Thanks

[ETHANol]

intersectBed from Bedtools. Can that do what you are trying to do?

[swNGS]

I'm having another stab at this question, and I've tried to visually represent it in an attempt to explain better what I am trying to do:

I have a library that is comprised of thousands of overlapping amplicons.
I have the genomic coordinates of all the amplicons that go to make the library.

I want to be able to count the number of occurrences of paired reads that are derived from the original amplicons. The rationale for this is to look for amplicon drop out in a single sample relative to other multiplexed samples as indicative of a possible variant under a site involved in generating the original amplicon (probe binding site or restriction endonuclease recognition site)

The following image attempts to explain it:
Green track is the expected amplicons.
Blue is the region of interest



Amplicon1 = chr19:11203020-11203404
Amplicon2 = chr19:11203059-11203231
Amplicon3 = chr19:11202984-11203057

I want to be able to count the aligned read pairs that represent each amplicon.
If no amplicons overlapped, I could probably do this relatively easily , but I'm slightly stumped about what to do with situations like amplicons 1& 2... So in that circumstance, I dont want reads from amplicon2 counting towards amplicon1.

The end result should be something like:

Genomic_coord AmpliconID Ocurrances
chr19:11203020-11203404 Amplicon1 153
chr19:11203059-11203231 Amplicon2 15
chr19:11202984-11203057 Amplicon3 48

While I can find plenty of tools to produce summary info on insert size ranges, I can find anything to specifically do this.

I think intersectBed is probably a place to start as suggested above, but cant get my head around the caveats that I listed. I was hoping there would be other helpful suggestions out there.

Thanks,

Chris

[frozenlyse]

In R, the GenomicRanges package has the function countOverlaps which allows you to specify type="within" - that way the query must be completely inside the subject.

[swNGS]

How would that work if I have two amplicons where one lies entirely within another?
I'm not sure if the image has uploaded correctly so here is a link

http://s14.postimage.org/rbmswg41r/igv_snapshot_Main.png

[frozenlyse]

Ah sorry I didn't notice one was inside the other.

One (hacky) thing you could do is to subtract the counts of amplicon 2 from amplicon 3. Another is that if the reads from amplicons will always start and end at defined positions, you could use type="equal", so that way only reads that start and end exactly at your amplicon start/ends will be counted.

[Chipper]

I would use a short perl thing with the fragment start and length as a hash key:

samtools view file.bam | perl -e ' while (<>) {split; $id="$_[2]:$_[3]:$_[8]";$hash{$id}++;}foreach (sort keys %hash){print $_ ."\t".$hash{$_}."\n";}' > fragments.count