Seqanswers Leaderboard Ad

Collapse
This topic is closed.
X
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • jullee
    Member
    • Apr 2014
    • 19

    How to compare Read IDs from different Bam files

    Hi

    I have aligned (and filtered) my paired-end sequence data to two alternative reference genomes (e.g. different species). I would like to get a list of reads that align to both references and a list of reads that only align to one. My approach so far is as follows:

    1) convert my final bam files to sam files

    samtools view -h -o Ref1.sam Ref1.bam

    samtools view -h -o Ref2.sam Ref2.bam

    2) extract column 1 (sequence IDs) from both files (additional grep step to remove the header)

    awk '{print $1}' Ref1.sam > Ref1Ids
    grep -v '^@' Ref1Ids > Ref1Ids_Final

    (repeat for Ref2 sam file)

    3) ?

    Now I'm thinking about using basic linux commands to compare the text files listing sequence IDs. But I'm having trouble figuring out how. I think the "comm" command requires that the files be "sorted" but when I use the "sort" command on my files, the resulting order of read IDs in the files doesn't make sense (e.g. read names are composed of both strings and numbers; the strings are constant but the numbers vary and are out of a logical order following the sort command).

    Anyone have any suggestions? Or a better way to get a list of read IDs that differ between two bam alignment files?

    Thanks!
  • lindenb
    Senior Member
    • Apr 2010
    • 143

    #2
    I wrote a tool to compare two BAMS: https://github.com/lindenb/jvarkit/wiki/CmpBams

    Comment

    • kmcarr
      Senior Member
      • May 2008
      • 1181

      #3
      Hi Jullee,

      Here is how I would do it:

      Code:
      samtools view -F 4 Ref1.bam | cut -f1 | sort -u > Ref1_reads.txt
      samtools view -F 4 Ref2.bam | cut -f1 | sort -u > Ref2_reads.txt
      
      comm -12 Ref1_reads.txt Ref2_reads.txt > AlignToBoth.txt
      comm -23 Ref1_reads.txt Ref2_reads.txt > AlignToRef1Only.txt
      comm -13 Ref1_reads.txt Ref2_reads.txt > AlignToRef2Only.txt
      Let me break this down:

      The samtools command is converting the the BAM to SAM format but filtering out any reads which are not aligned to the reference (-F 4).

      You should not output the header information (-h) since it is meaningless for your purposes (and you got rid of it later in your pipeline anyway).

      I do not bother saving the intermediate .sam file since I really don't need to, just pipe the output into the next step.

      Just take the read ID in column 1 with cut. (You used awk for this but whatever you like.)

      Sort the ID's and only keep one, unique copy of each (-u). Save the output.

      You were concerned about the order output by sort; don't be. sort in this case is outputting the read IDs in lexical (dictionary) order. This is the order which is required by comm.

      The comm commands are pretty straightforward, outputting just one of the three columns in each run. You could run comm just once and then separate the columns but I find this method a little easier.

      ONE BIG NOTE: This pipeline will collapse the paired read IDs into a single entry meaning that only one of the reads from a pair need align to either of the references to be counted.

      Comment

      • jullee
        Member
        • Apr 2014
        • 19

        #4
        Thank you kmcarr and lindenb for the helpful replies! I will try these out.

        Comment

        Latest Articles

        Collapse

        • seqadmin
          New Genomics Tools and Methods Shared at AGBT 2025
          by seqadmin


          This year’s Advances in Genome Biology and Technology (AGBT) General Meeting commemorated the 25th anniversary of the event at its original venue on Marco Island, Florida. While this year’s event didn’t include high-profile musical performances, the industry announcements and cutting-edge research still drew the attention of leading scientists.

          The Headliner
          The biggest announcement was Roche stepping back into the sequencing platform market. In the years since...
          03-03-2025, 01:39 PM
        • seqadmin
          Investigating the Gut Microbiome Through Diet and Spatial Biology
          by seqadmin




          The human gut contains trillions of microorganisms that impact digestion, immune functions, and overall health1. Despite major breakthroughs, we’re only beginning to understand the full extent of the microbiome’s influence on health and disease. Advances in next-generation sequencing and spatial biology have opened new windows into this complex environment, yet many questions remain. This article highlights two recent studies exploring how diet influences microbial...
          02-24-2025, 06:31 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, Today, 05:03 AM
        0 responses
        15 views
        0 reactions
        Last Post seqadmin  
        Started by seqadmin, Yesterday, 07:27 AM
        0 responses
        12 views
        0 reactions
        Last Post seqadmin  
        Started by seqadmin, 03-18-2025, 12:50 PM
        0 responses
        14 views
        0 reactions
        Last Post seqadmin  
        Started by seqadmin, 03-03-2025, 01:15 PM
        0 responses
        185 views
        0 reactions
        Last Post seqadmin  
        Working...