-2
\$\begingroup\$

I was able to make a program that parse my samtools output file into a table that can be viewed easily in excel. While I was able to get the final result, is there a recommended improvements that be done for my code. I am trying to improve my python skills and transferring over from C and C++.

The strategy I used was to split the data I need based on the "+" sign and make it into an array. Then select the elements of the array that were the information that I needed and write it to my file.

My example input:

15051874 + 0 in total (QC-passed reads + QC-failed reads)
1998052 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
13457366 + 0 mapped (89.41% : N/A)
13053822 + 0 paired in sequencing
6526911 + 0 read1
6526911 + 0 read2
10670914 + 0 properly paired (81.75% : N/A)
10947288 + 0 with itself and mate mapped
512026 + 0 singletons (3.92% : N/A)
41524 + 0 with mate mapped to a different chr
31302 + 0 with mate mapped to a different chr (mapQ>=5)

My output:

FileName    Total   Secondary   Supplementary   duplicates  mapped  paired in sequencing    read1   read2   properly paired with itself and mate mapped singletons  with mate mapped to a different chr with mate mapped to a different chr (mapQ>=5)
10_HK_S22.merged.samtools.flag.txt  26541257    2332283 0   0   22895440    24208974    12104487    12104487    19003826    19632880    930277  69030   52261

My Program:

outFile = open("output.count.txt", "w+")

#windows platform add the r
os.chdir(r"Susceptible\featurecounts")


#open the output file to be able to write output.
outFile.write("FileName\tTotal\tSecondary\tSupplementary\tduplicates\tmapped\tpaired in sequencing\tread1\t"
    "read2\tproperly paired\twith itself and mate mapped\tsingletons\twith mate mapped to a different chr\twith mate mapped to a different chr (mapQ>=5)\n")


#Iterate through files in directory with the following ending
for file in glob.glob(".flag.txt"):
    #open file after retrieving the name.
    with open(file, 'r') as counts_file:
        #empty list/array for storing the outputs
        list = []
        #add the file name to array. 
        list.append(file)
        #get values from output file.
        for line in counts_file:
            list.append(line.split('+')[0])

        #write list to file
        for item in list:
            outFile.write("%s\t" % item)

        #write a newline
        outFile.write("\n")


#close the output file      
outFile.close()
\$\endgroup\$
  • 1
    \$\begingroup\$ Your example output doesn't seem to correspond to your example input? \$\endgroup\$ – 200_success Oct 2 '18 at 18:17
  • \$\begingroup\$ Please take a look at the How do I ask a good question? page. \$\endgroup\$ – code_dredd Oct 2 '18 at 18:31
  • \$\begingroup\$ @200_success would have been better to provide a google link to the example? \$\endgroup\$ – user1238097 Oct 2 '18 at 18:48
  • \$\begingroup\$ No, I just expect your numbers to be consistent between your example input and example output. \$\endgroup\$ – 200_success Oct 2 '18 at 18:51
  • \$\begingroup\$ @200_success Ah yeh, I see what you mean completely now. \$\endgroup\$ – user1238097 Oct 2 '18 at 18:52
0
\$\begingroup\$

Use with ... as ...: statements to open files, and automatically close them. Then you don't have to clutter up your program with explicit close statements.

outFile = open("output.count.txt", "w+")

# ... code here

#close the output file      
outFile.close()   

Becomes:

with open("output.count.txt", "w+") as outFile:

    # ... code here

This is ugly and unreadable:

outFile.write("FileName\tTotal\tSecondary\tSupplementary\tduplicates\tmapped\tpaired in sequencing\tread1\t"
    "read2\tproperly paired\twith itself and mate mapped\tsingletons\twith mate mapped to a different chr\twith mate mapped to a different chr (mapQ>=5)\n")

The \t runs into the next field name, so the eye sees "tTotal". It would be better to actually list your field names in a readable form, and let the computer properly separate them:

fields = ["FileName", "Total", "Secondary", "Supplementary", "duplicates", "mapped",
          "paired in sequencing", "read1", "read2", "properly paired",
          "with itself and mate mapped", "singletons", "with mate mapped to a different chr",
          "with mate mapped to a different chr (mapQ>=5)"]

outFile.write("\t".join(fields) + '\n')

Looping through one iterable, processing each one and creating a new list be often done cleaner using list comprehension:

    list = []
    #add the file name to array. 
    list.append(file)
    #get values from output file.
    for line in counts_file:
        list.append(line.split('+')[0])

Could become (without the "file" at the start of the list):

    values = [ line.split('+')[0] for line in counts_file ]

But you take the resulting list and add a \t character between each value, so maybe instead:

    values = "\t".join( line.split('+')[0] for line in counts_file )

Now, you want to print out the values to the outFile, with the file at the start. f-strings are a new feature in Python. They let you format a string with local variables interpolated into the string. This makes it easy:

    outFile.write(f"{file}\t{values}\n")

As a bonus, each line doesn't end with a trailing tab character.


Resulting code would be something like:

with open("output.count.txt", "w+") as outFile:

    fields = ["FileName", "Total", "Secondary", "Supplementary", "duplicates", "mapped",
              "paired in sequencing", "read1", "read2", "properly paired",
              "with itself and mate mapped", "singletons", "with mate mapped to a different chr",
              "with mate mapped to a different chr (mapQ>=5)"]

    outFile.write("\t".join(fields) + '\n')

    for file in glob.glob(".flag.txt"):
        with open(file, 'r') as counts_file:
           values = "\t".join( line.split('+')[0] for line in counts_file )
           outFile.write(f"{file}\t{values}\n")
\$\endgroup\$
  • 1
    \$\begingroup\$ Note: string interpolation is present in python from version 3.6 \$\endgroup\$ – MaLiN2223 Oct 3 '18 at 8:40

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