Script to find and print common IDs between two files working but it's not optimal

Gf.Ena 0
Viewed 44

I have a code that is working and does what I want, but it is extremely slow. It takes 1 or 2 days depending on the size of the input files. I know that there are alternatives that can be almost instant and that my code is slow because it's a recursive grep. I wrote another code in python that works as intended and is almost instant, but it does not print everything I need. What I need is the common IDs between two files, and I want it to print the whole line. My python script does not do that, while the bash does it but it's too much slow.

This is my code in bash:

awk '{print $2}'  file1.bim > sites.txt 

for snp in `cat sites.txt`
do

grep -w $snp file2.bim >> file1_2_shared.txt

done

This is my code in python:

#!/usr/bin/env python3
import sys
argv1=sys.argv[1]   #argv1 is the first .bim file
argv2=sys.argv[2]   #argv2 is the second .bim file
argv3=sys.argv[3]   #argv3 is the output .txt file name

def printcommonSNPs(inputbim1,inputbim2,outputtxt):
    bim1 = open(inputbim1, "r")
    bim2 = open(inputbim2, "r")
    output = open(outputtxt,"w")
    
    snps1 = []
    line1 = bim1.readline()
    line1 = line1.split()
    snps1.append(line1[1])
    for line1 in bim1:
        line1 = line1.split()
        snps1.append(line1[1])
    bim1.close()
        
    snps2 = []
    line2 = bim2.readline()
    line2 = line2.split()
    snps2.append(line2[1])
    for line2 in bim2:
        line2 = line2.split()
        snps2.append(line2[1])
    bim2.close()
    common=[]
    common = list(set(snps1).intersection(snps2))
    
    for SNP in common:
        print(SNP, file=output)
    
printcommonSNPs(argv1,argv2,argv3)

My .bim input files are made this way:

1   1:891021    0   891021  G   A
1   1:903426    0   903426  T   C
1   1:949654    0   949654  A   G

I would appreciate suggestions on what I could do to make it quick in bash (I suspect I can use an awk script, but I tried awk 'FNR==NR {map[$2]=$2; next} {print $2, map[$2]}' file1.bim file2.bim > Roma_sets_shared_sites.txt and it simply prints every line, so it's not working as I need), or how could I tell to print the whole line in python3.

grep awk bash for-loop
1 Answers

It looks as if the problem can be solved like this:

grep -w -f <(awk '{ print $2 }' file1.bim) file2.bim

The identifiers (field $2) from file1.bim are to be treated as patterns to grep for in file2.bim. GNU grep takes a -f file argument which gives a list of patterns, one per line. We use <() process substitution in place of a file. It looks as if the -w option individually applies to the -f patterns.

This won't have the same output as your shell script if there are duplicate IDs in file1.bim. If the same pattern occurs more than once, that's the same as one instance. And of course the order is different. Grepping the entire second file for one identifier and hen the next and next, produces the matches in a different order. If that order has to be reproduced, it will take extra work.

usermpf50
Related