Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • coverageBed counts not the same as input bed

    When I use a sorted bed file with unique 56,548 targets in it and process it in bedtools 2.25.0 with the below command:

    Code:
    coverageBed -d -g /home/cmccabe/Desktop/NGS/bedtools2-25.0/genomes/human.hg19.Ion.genome.txt -a /home/cmccabe/Desktop/NGS/bed/bedtools/bedtools_unix_sorted_xgen_targets.bed -b /home/cmccabe/Desktop/NGS/IonXpress_004.R_2015_10_29_13_11_15_user_Proton-32-Lurie_Inh_Disease_151029_newheader.bam > /home/cmccabe/Desktop/NGS/IonXpress_004.R_2015_10_29_13_11_15_user_Proton-32-Lurie_Inh_Disease_151029_newheader_basecounts.bed
    I get 56,543 unique targets returned. The sort order is lexical in the bed and human in the bam so I use the -g option to specific the ordering of the bam. The 56,543 was confirmed in perl and awk. Any ideas why the two #'s are not the same? Thank you .

  • #2
    Perhaps those 5 targets are not covered by reads in your BAM? Did you look at the BAM in IGV?

    Comment


    • #3
      I will try that, but shouldn't bedtools output that region with all zero's.

      There are other regions in the coverageBed output where the regions is zero.

      For example each one of my targets is 120 bases and the coverage per base is calculated. So if there was no coverage in those 120 bases, there is no output in the file? Thank you .

      Comment


      • #4
        I see.

        Just for the heck of it can you do a "sort/uniq" on the names to see if those 5 are for some reason repeated in the file?

        Edit: Just thought that the intervals may be repeated (instead of just names) as another possibility.
        Last edited by GenoMax; 11-12-2015, 09:24 AM.

        Comment


        • #5
          If I do
          Code:
          wc -l
          on the original bed, I get 56,548, but if I do:
          Code:
          awk '$4!=d{c++;d=$4}END{print c}' file
          or perl -nae '$seen{$F[3]}++;
          END{
          print "There are ", scalar keys %seen, " unique fourth fields\n";
          }' originalbed [/CODE]

          I get 56,543

          However, in excel there are 56,548. Since the file was sent to me I think I am just going to re-do it without using excel (which is very problematic). Thank you .

          Comment


          • #6
            5 blank lines at the end of that file?

            Comment


            • #7
              I have put the bed file on box.net as it is too large to attach:



              When I do a wc-l test.bed I get 56,548

              However, when I use awk or perl I get 56,543 and that is also what bedtools outputs.

              Code:
              cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ wc -l test.bed
              56548 test.bed
              cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ awk '$4!=d{c++;d=$4}END{print c}' test.bed
              56543
              cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed/bedtools$ perl -nae '$seen{$F[3]}++;
              >     END{
              >         print "There are ", scalar keys %seen, " unique fourth fields\n";
              >     }' test.bed
              There are 56543 unique fourth fields
              I can not seem to figure out why there are 5 missing. Thank you .

              edit:
              you were right @Genomax they were non-unique and excel didn't mark them:
              Code:
              4748d4747
              < chr11	47270217	47270425	chr11:47270217-47270425	unknown-1062|gc=64.9
              4970d4968
              < chr11	5248271	5248449	chr11:5248271-5248449	HBB-283|gc=55.1
              24883d24880
              < chr19	13010118	13010237	chr19:13010118-13010237	SYCE2-864|gc=47.9
              33027d33023
              < chr22	38153605	38154160	chr22:38153605-38154160	TRIOBP-610|gc=68.6
              54957d54952
              < chrX	33357316	33359011	chrX:33357316-33359011	DMD-581|gc=33.7
              but fgrep did
              Code:
              $ fgrep 'chr11:47270217-47270425' test.bed
              chr11	47270217	47270425	chr11:47270217-47270425	ACP2-1062|gc=64.9
              chr11	47270217	47270425	chr11:47270217-47270425	unknown-1062|gc=64.9
              Thank you
              Last edited by cmccabe; 11-12-2015, 04:04 PM. Reason: added edit

              Comment


              • #8
                Ah. It is good to have a logical explanation.

                Comment

                Latest Articles

                Collapse

                • 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
                • seqadmin
                  Choosing Between NGS and qPCR
                  by seqadmin



                  Next-generation sequencing (NGS) and quantitative polymerase chain reaction (qPCR) are essential techniques for investigating the genome, transcriptome, and epigenome. In many cases, choosing the appropriate technique is straightforward, but in others, it can be more challenging to determine the most effective option. A simple distinction is that smaller, more focused projects are typically better suited for qPCR, while larger, more complex datasets benefit from NGS. However,...
                  10-18-2024, 07:11 AM

                ad_right_rmr

                Collapse

                News

                Collapse

                Topics Statistics Last Post
                Started by seqadmin, 11-08-2024, 11:09 AM
                0 responses
                223 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 11-08-2024, 06:13 AM
                0 responses
                165 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 11-01-2024, 06:09 AM
                0 responses
                80 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 10-30-2024, 05:31 AM
                0 responses
                27 views
                0 likes
                Last Post seqadmin  
                Working...
                X