Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

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

    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

  • #2
    intersectBed from Bedtools. Can that do what you are trying to do?
    --------------
    Ethan

    Comment


    • #3
      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

      Comment


      • #4
        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.

        Comment


        • #5
          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/i...pshot_Main.png

          Comment


          • #6
            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.

            Comment


            • #7
              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

              Comment

              Latest Articles

              Collapse

              • seqadmin
                Latest Developments in Precision Medicine
                by seqadmin



                Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

                Somatic Genomics
                “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
                Yesterday, 01:16 PM
              • seqadmin
                Recent Advances in Sequencing Analysis Tools
                by seqadmin


                The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
                05-06-2024, 07:48 AM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, Yesterday, 07:15 AM
              0 responses
              13 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 05-23-2024, 10:28 AM
              0 responses
              17 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 05-23-2024, 07:35 AM
              0 responses
              20 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 05-22-2024, 02:06 PM
              0 responses
              10 views
              0 likes
              Last Post seqadmin  
              Working...
              X