0
\$\begingroup\$

First of all, please note that the input (.vcf) file is very large (>60 GB).

So I am reading the .vcf files line by line, but since I have to check for duplicate of two columns present in the input vcf file, I extract those two columns from the input file to a text file ("dup_posID.txt") and then check for duplicates, and remove the duplicates and save it to a new text file ("dup_posID_re.txt").
Then I create a new vcf file (.redp.vcf) by comparing the recoded text file.

The code works fine and am getting the desired result. But its taking a lot of time. more that 24 hours. How can I reduce that time?

vcf_File_name = "chr7_dup_fin"

with  open(vcf_File_name+".vcf","r") as f, open(vcf_File_name+"dup_posID.txt", "w+") as D:

    for v_line in  f:

        if v_line[0] != '#':
            v_cols = re.split('[ ]+|[   ]+',v_line.strip())
            r= v_cols[2].strip() +"\t"+ v_cols[1].strip() + "\n"
            req1 += r
            #print(req1)
    D.write(req1)

D1=open(vcf_File_name+"dup_posID.txt","r")

D2 = open (vcf_File_name+"dup_posID_re.txt", "w+")

lines1 = D1.readlines()

for line1 in lines1:
    for line1 in lines1:
        cols= re.split('[   ]', line1.strip())
        try:
            dup[cols[0]] += 1 
            #print cols[0]
        except:

            try:
                dup_pos[cols[1]] += 1
                #print "Duplicate Position: "+cols[1]
            except:
                dup_pos[cols[1]] = 0
                dup[cols[0]] = 0


                p = cols[0].strip()+"\t"+cols[1].strip()+"\n"
                q += p
D2.write(q)

D2.close()




with  open(vcf_File_name+".vcf","r") as f, open(vcf_File_name+".redp.vcf","w+") as f1 :
    D2=open(vcf_File_name+"dup_posID_re.txt", "r")
    lines2=D2.readlines()
    for v_line in  f:

        if v_line[0] != '#':
            v_cols = re.split('[ ]+|[   ]+',v_line.strip())

            st = v_cols[2]+"\t"+v_cols[1]+"\n"
            if st in lines2:
                lines2.remove(st)
                s = "\t".join(v_cols)
                s += "\n"

                #req += s
                f1.write(s)
                #f1.write("\n")
        else:
            s = v_line
            f1.write(s)
    f1.close()

The input .vcf file has 2512 columns and is tab separated. Example:

2   26813   rs574377045  C  A   100 PASS    AC=2;AF=0.000399361;AN=5008;NS=2504;DP=18198;EAS_AF=0.001;AMR_AF=0.0014;AFR_AF=0;EUR_AF=0;SAS_AF=0;VT=SNP;AA=C GT   0|0
2   26816   rs535194004  C  T   100 PASS   AC=1;AF=0.000199681;AN=5008;NS=2504;DP=18346;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;EUR_AF=0;SAS_AF=0;VT=SNP;AA=C      GT   0|0

Here, one line extends till 2512 columns.

\$\endgroup\$
  • 1
    \$\begingroup\$ Please include all the code you're using. req1 is missing the first assignment for example. If you could give a small description of what your code is performing then it'd help a lot too. \$\endgroup\$ – Peilonrayz Apr 24 '18 at 10:02
  • 1
    \$\begingroup\$ Welcome to Code Review! This question is incomplete. To help reviewers give you better answers, please add sufficient context to your question. The more you tell us about what your code does and what the purpose of doing that is, the easier it will be for reviewers to help you. The current title states your concerns about the code; it needs an edit to simply state the task; see How to get the best value out of Code Review: Asking Questions for guidance on writing good question titles. \$\endgroup\$ – Toby Speight Apr 24 '18 at 11:52
2
\$\begingroup\$
  1. Stop reading and writing to files so much. You seem to be able to read the entire file into RAM since D1.readlines() Doesn't cause your program to explode.
  2. Create a function to read the file, split using your regex and return the columns.
  3. Rather than using a regex I'd use str.split().

This can create:

def read_file(file_name):
    with open(file_name) as f:
        for line in  f:
            if line.startswith('#'):
                yield tuple(), line
            else:
                columns = line.strip().split()
                yield tuple(columns[1:3]), line
  1. You're removing duplicate ID's with your try excepts. You can do this with a set.
  2. You can make this work with the above by filtering empty columns, as they're commented out code.

Creating:

def all_columns(file_name):
    return set(
        cols
        for cols, _ in read_file(file_name)
        if cols
    )
  1. Now we can get all the column names without duplicates by doing columns = all_columns(file_name + '.vcf').
  2. Open the file to write to.
  3. Read the original file again, so that we can filter the columns.
  4. Write all lines that are commented out.
  5. If the column is in all the collumns, then we will:

    1. Write the line to the output.
    2. Remove the column from the set.

Creating:

if __name__ == '__main__':
    file_name = 'chr7_dup_fin'

    columns = all_columns(file_name + '.vcf')
    with open(file_name + '.redp.vcf', 'w+') as output:
        for cols, line in read_file(file_name + '.vcf'):
            if not cols:
                output.write(line)
            elif cols in columns:
                output.write(line)
                columns.remove(cols)
\$\endgroup\$
2
\$\begingroup\$

Some notes:

  1. Run it through pycodestyle to learn how your style differs (at least in machine detectable ways) from idiomatic Python. This makes code not only easier to review, but also easier for seasoned Python developers to read and maintain. I would recommend not silencing any of the complaints, and instead fixing every one of them. It quickly becomes second nature to write idiomatic Python that way.
  2. Abbr. var nams hindr rdability, makin it hardr to rd. Unless D1 is absolutely clear from context to anyone in your field I would rename it until it is. The same of course goes for f, which could be semi-tolerable if your script only ever processed one file, but is baffling if you're processing several. "More typing" is not an argument against this, because there are editors which will happily auto-complete any identifier for you.
  3. You have a nested loop where both loops use the same instance variable. That means that the outer loop variable is useless.
  4. You have huge files you want to process, so you really, really don't want to read whole files before starting processing. The way to avoid that is to process each line as you read it, as in for line in file_pointer rather than lines = file.readlines() + a loop.
  5. The main parts of this should be split into functions, called by a short main function which is triggered only if __name__ == '__main__': - this rather unfortunate pattern is idiomatic Python, and allows reuse of the functions and objects in your module.
  6. You have completely open ended except clauses, which is almost universally a bad idea. Exception handling is expensive, and you'll be much better off checking for the existence of an index than throwing an exception and catching it again immediately. In general, the anti-pattern is to use exceptions for flow control.
  7. There are several patterns in use which may be fine for most processing, but which are not optimal for large amounts of data, such as splitting and joining of strings, regular expressions, and accumulating strings and then writing once rather than just writing every substring. From a cursory look it seems that your processing would be very much quicker if done with a database, if at all possible.
  8. You are taking care to close some files by using with open(), but the second D2 is not closed.
  9. Reusing variables for different things also makes code really hard to read, and it a common misconception that it is more efficient. Maybe in some niche corner cases, but for this kind of processing you will be much better off with evidence-based optimisation - prove where the code is slow with a profiler and fix those parts.
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.