This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How Do I Do The Intersection Of Two Sam Files

I have two sam files produced by alignment to two different versions of the reference genome. Is there a way to get the intersection of the sam files so I get a sam with the sequences that have been mapped to both? (So, I would ideally want to discard any sequences that are unmapped in an alignment to at least one of the genomes.)

sam
• 4,345 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by KCC KCC

2 answers

You can try using this python script. I am not 100% sure it will work though. I am at home right now and don't have access to any of my SAM files to test it out:

import sys

inFileA = open(sys.argv[1],'r')
inFileB = open(sys.argv[2],'r')

reads = {}
for line in inFileA:
    if line[0] != "@": 
        data = line.strip().split('\t')
        name = data[0]
        flag = int(data[1])
        if not flag & 4:
            reads[name] = True

for line in inFileB:
    if line[0] != "@": 
        data = line.strip().split('\t')
        name = data[0]
        flag = int(data[1])
        if not flag & 4:
            if reads.has_key(name):
                print name

Use it by: python theScript.py fileA.sam fileB.sam > listOfReadNames

• 0 views • link
written by Damian Kao Damian Kao
  • use samtools view with option '-f' to dump the sequences mapped on both assemblies
  • pipe in cut to extract the names of the reads
  • pipe in sort to order the names
  • redirect to file

use comm or join to list the common names in both files

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this