Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Sort/create list of locations with highest # of reads?

    Hello,

    Hopefully there is a simple way to do this that I am unaware of.

    I have fastq and bam files (100bp paired end Illumina). In the alignment files there are lots of areas with one or a few reads spread throughout the genome.

    What I would like to do is either sort results to visualize, for example in Integrative Genomics Viewer, or generate a list that will show in descending order the locations that have the greatest number of reads or the greatest coverage. This would most likely be within a window, say 50bp.

    For example:

    chr21:33,031,597 has 58 reads
    chr21:33,031,200 has 56 reads
    chr21:33,025,000 has 42 reads

    and so on.

    Any advice or guidance on different ways I might achieve this is very appreciated.

    Thank you

  • #2
    Check out samtools ...

    $ samtools depth

    Usage: samtools depth [options] in1.bam [in2.bam [...]]
    Options:
    -b <bed> list of positions or regions
    -f <list> list of input BAM filenames, one per line [null]
    -l <int> minQLen
    -q <int> base quality threshold
    -Q <int> mapping quality threshold
    -r <chr:from-to> region

    Example:
    $ samtools depth -r chr1:1-100000 1221792.bam | head
    chr1 10033 1
    chr1 10034 1
    chr1 10035 1
    chr1 10036 1
    chr1 10037 1
    chr1 10038 1
    chr1 10039 1
    chr1 10040 1
    chr1 10041 1
    chr1 10042 1


    Use "awk" to filter on column 3 and "sort" on column 3 and "tail" it.
    Like this ...
    samtools depth 1221792.bam | awk '{if ($3>300) print $0}' | sort -k3n | tail
    Last edited by Richard Finney; 06-06-2014, 06:43 AM.

    Comment


    • #3
      Thanks so much Richard,

      So if I use:

      samtools depth 1221792.bam | awk '{if ($3>50) print $0}' | sort -k3n | tail

      for example, this will give me a list of all locations with more than 50 reads, correct?

      But no 'window,' or in other words, it will likely print the individual locations for each block of reads, rather than one result for each bin or chunk of 50 bases or something like that. Meaning, at chr1 10038 +/- 25bp, there are max 50 reads.

      Certainly a start, thank you for the help.

      Comment


      • #4
        Code:
        samtools view -F 4 sorted.bam | cut -f 3,4 | uniq -c | sort -nr | head -n 20 > most_common.txt
        That'll give you the top 20 most common start sites for your reads.

        Comment

        Latest Articles

        Collapse

        • seqadmin
          The Impact of AI in Genomic Medicine
          by seqadmin



          Article Coming Soon......
          Today, 02:07 PM
        • seqadmin
          Multiomics Techniques Advancing Disease Research
          by seqadmin


          New and advanced multiomics tools and technologies have opened new avenues of research and markedly enhanced various disciplines such as disease research and precision medicine1. The practice of merging diverse data from various ‘omes increasingly provides a more holistic understanding of biological systems. As Maddison Masaeli, Co-Founder and CEO at Deepcell, aptly noted, “You can't explain biology in its complex form with one modality.”

          A major leap in the field has
          ...
          02-08-2024, 06:33 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, 02-23-2024, 04:11 PM
        0 responses
        31 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 02-21-2024, 08:52 AM
        0 responses
        46 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 02-20-2024, 08:57 AM
        0 responses
        36 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 02-14-2024, 09:19 AM
        0 responses
        63 views
        0 likes
        Last Post seqadmin  
        Working...
        X