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
                Recent Advances in Sequencing Technologies
                by seqadmin







                Innovations in next-generation sequencing technologies and techniques are driving more precise and comprehensive exploration of complex biological systems. Current advancements include improved accessibility for long-read sequencing and significant progress in single-cell and 3D genomics. This article explores some of the most impactful developments in the field over the past year.

                Long-Read Sequencing
                Long-read sequencing has...
                12-02-2024, 01:49 PM
              • seqadmin
                Genetic Variation in Immunogenetics and Antibody Diversity
                by seqadmin



                The field of immunogenetics explores how genetic variations influence immune responses and susceptibility to disease. In a recent SEQanswers webinar, Oscar Rodriguez, Ph.D., Postdoctoral Researcher at the University of Louisville, and Ruben Martínez Barricarte, Ph.D., Assistant Professor of Medicine at Vanderbilt University, shared recent advancements in immunogenetics. This article discusses their research on genetic variation in antibody loci, antibody production processes,...
                11-06-2024, 07:24 PM

              ad_right_rmr

              Collapse

              News

              Collapse

              Topics Statistics Last Post
              Started by seqadmin, 12-02-2024, 09:29 AM
              0 responses
              145 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 12-02-2024, 09:06 AM
              0 responses
              51 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 12-02-2024, 08:03 AM
              0 responses
              42 views
              0 likes
              Last Post seqadmin  
              Started by seqadmin, 11-22-2024, 07:36 AM
              0 responses
              72 views
              0 likes
              Last Post seqadmin  
              Working...
              X