Announcement

Collapse
No announcement yet.

extracting identical reads from bam file

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • extracting identical reads from bam file

    Hello,
    I have an alignment file of smallRNAs in bam format.
    How can I extact from this file all the identical reads in a single sequence but mantaining the reads count?
    Thank you for the help!

  • #2
    Are you trying to match a specific sequence ?

    samtools view x.bam | grep TTTTCTGCCTGTTGGGCTGGAG | awk '{print $10}' | uniq -c


    All reads with same sequence as another read, try this ...

    samtools view x.bam | awk '{print $10}' | sort --buffer-size=20G | uniq -c | awk '{if ($1!=1) print $0}'


    where x.bam is you bam file.
    fine tune the --buffer-size parameter to sort

    Comment


    • #3
      It would take longer, but swap the second awk in Richard's command line with

      Code:
      sort -nr
      To get the most common reads in order, starting with the most abundant. Tack a

      Code:
      | head -n 100
      to get only the top 100

      Comment


      • #4
        thank you! it worked!
        can I do the same from a fastq file?

        Comment


        • #5
          Yes, you might want to use awk to only print every 4th line.
          Note 1) the "mod" operator and 2) "line count" intrinsic variable in awk.
          Perl/python/C/java if you prefer can address the issue of filtering for only the sequence also.

          Comment


          • #6
            After I extracted the identical reads in a single sequence from the bam file I aligned them again to the genome. When I use igv to visualize the alignment now all the sequences are mapping in sense orientation.. even those sequences that are supposed to be antisense to the genome are shown in sense orientation. Why is that?

            Comment


            • #7
              The sequence for the read in a bam files may be reverse complemented to align to the reference.
              Reads are supposed to be properly noted as reversed in the bitwise flags field in a line/entry of sam/bam.

              You may wish to interrogate this flag for special processing.

              Comment

              Working...
              X