Command-line Arguments

If you have not already, I encourage you to copy the "new_py.py" script into your $PATH and then execute it with the -a argument to start a new script with argparse:

$ ./new_py.py -a test
Done, see new script "test.py."

If you check out the new script, it has a get_args function that will show you how to create named arguments for strings, integers, booleans, and positional arguments:

$ cat -n test.py
     1    #!/usr/bin/env python3
     2
     3    import argparse
     4    import os
     5    import sys
     6
     7    # --------------------------------------------------
     8    def get_args():
     9        parser = argparse.ArgumentParser(description='Argparse Python script')
    10        parser.add_argument('positional', metavar='str', help='A positional argument')
    11        parser.add_argument('-a', '--arg', help='A named string argument',
    12                            metavar='str', type=str, default='')
    13        parser.add_argument('-i', '--int', help='A named integer argument',
    14                            metavar='int', type=int, default=0)
    15        parser.add_argument('-f', '--flag', help='A boolean flag',
    16                            action='store_true')
    17        return parser.parse_args()
    18
    19    # --------------------------------------------------
    20    def main():
    21        args = get_args()
    22        str_arg = args.arg
    23        int_arg = args.int
    24        flag_arg = args.flag
    25        pos_arg = args.positional
    26
    27        print('str_arg = "{}"'.format(str_arg))
    28        print('int_arg = "{}"'.format(int_arg))
    29        print('flag_arg = "{}"'.format(flag_arg))
    30        print('positional = "{}"'.format(pos_arg))
    31
    32    # --------------------------------------------------
    33    if __name__ == '__main__':
    34        main()

If you run without any arguments, you get a nice usage statement:

$ ./test.py
usage: test.py [-h] [-a str] [-i int] [-f] str
test.py: error: the following arguments are required: str
$ ./test.py foobar
str_arg = ""
int_arg = "0"
flag_arg = "False"
positional = "foobar"
$ ./test.py -a XYZ -i 42 -f foobar
str_arg = "XYZ"
int_arg = "42"
flag_arg = "True"

Please RTFM (https://docs.python.org/3/library/argparse.html\ to find all the other Fine things you can do with argparse.

CSV Files

Delimited text files are a standard way to distribute non/semi-hierarchical data -- e.g., records that can be represented each on one line. (When you get into data that have relationships, e.g., parents/children, then structures like XML and JSON are more appropriate, which is not to say that people haven't sorely abused this venerable format, e.g., GFF3.) Let's first take a look at the csv module in Python (https://docs.python.org/3/library/csv.html\ to parse the output from Centrifuge (http://www.ccb.jhu.edu/software/centrifuge/.

For this, we'll use some data from a study from Yellowstone National Park (https://www.imicrobe.us/sample/view/1378\ a file ("YELLOWSTONE_SMPL_20723.sum") showing the taxonomy ID for each read it was able to classify and 2) a file ("YELLOWSTONE_SMPL_20723.tsv") of the complete taxonomy information for each taxonomy ID. One record from the first looks like this:

      readID: Yellowstone_READ_00007510
       seqID: cid|321327
       taxID: 321327
       score: 640000
2ndBestScore: 0
   hitLength: 815
 queryLength: 839
  numMatches: 1

One from the second looks like this:

          name: synthetic construct
         taxID: 32630
       taxRank: species
    genomeSize: 26537524
      numReads: 19
numUniqueReads: 19
     abundance: 0.0

Read Count by taxID

Let's write a program that shows a table of the number of records for each "taxID":

$ cat -n read_count_by_taxid.py
     1    #!/usr/bin/env python3
     2    """Counts by taxID"""
     3
     4    import csv
     5    import os
     6    import sys
     7    from collections import defaultdict
     8
     9    args = sys.argv[1:]
    10
    11    if len(args) != 1:
    12        print('Usage: {} SAMPLE.SUM'.format(os.path.basename(sys.argv[0])))
    13        sys.exit(1)
    14
    15    sum_file = args[0]
    16
    17    _, ext = os.path.splitext(sum_file)
    18    if not ext == '.sum':
    19        print('File extention "{}" is not ".sum"'.format(ext))
    20        sys.exit(1)
    21
    22    counts = defaultdict(int)
    23    with open(sum_file) as csvfile:
    24        reader = csv.DictReader(csvfile, delimiter='\t')
    25        for row in reader:
    26            taxID = row['taxID']
    27            counts[taxID] += 1
    28
    29    print('\t'.join(['count', 'taxID']))
    30    for taxID, count in counts.items():
    31        print('\t'.join([str(count), taxID]))

As always, it prints a "usage" statement when run with no arguments. It also uses the os.path.splitext function to get the file extension and make sure that it is ".sum." Finally, if the input looks OK, then it uses the csv.DictReader module to parse each record of the file into a dictionary:

$ ./read_count_by_taxid.py
Usage: read_count_by_taxid.py SAMPLE.SUM
$ ./read_count_by_taxid.py YELLOWSTONE_SMPL_20723.tsv
File extention ".tsv" is not ".sum"
$ ./read_count_by_taxid.py YELLOWSTONE_SMPL_20723.sum
count    taxID
80    321332
6432    321327
19    32630

That's a start, but most people would rather see the a species name rather than the NCBI taxonomy ID, so we'll need to go look up the taxIDs in the ".tsv" file:

$ cat -n read_count_by_tax_name.py
     1    #!/usr/bin/env python3
     2    """Counts by tax name"""
     3
     4    import csv
     5    import os
     6    import sys
     7    from collections import defaultdict
     8
     9    args = sys.argv[1:]
    10
    11    if len(args) != 1:
    12        print('Usage: {} SAMPLE.SUM'.format(os.path.basename(sys.argv[0])))
    13        sys.exit(1)
    14
    15    sum_file = args[0]
    16
    17    basename, ext = os.path.splitext(sum_file)
    18    if not ext == '.sum':
    19        print('File extention "{}" is not ".sum"'.format(ext))
    20        sys.exit(1)
    21
    22    tsv_file = basename + '.tsv'
    23    if not os.path.isfile(tsv_file):
    24        print('Cannot find expected TSV "{}"'.format(tsv_file))
    25        sys.exit(1)
    26
    27    tax_name = {}
    28    with open(tsv_file) as csvfile:
    29        reader = csv.DictReader(csvfile, delimiter='\t')
    30        for row in reader:
    31            tax_name[row['taxID']] = row['name']
    32
    33    counts = defaultdict(int)
    34    with open(sum_file) as csvfile:
    35        reader = csv.DictReader(csvfile, delimiter='\t')
    36        for row in reader:
    37            taxID = row['taxID']
    38            counts[taxID] += 1
    39
    40    print('\t'.join(['count', 'taxID']))
    41    for taxID, count in counts.items():
    42        name = tax_name.get(taxID) or 'NA'
    43        print('\t'.join([str(count), name]))
$ ./read_count_by_tax_name.py YELLOWSTONE_SMPL_20723.sum
count    taxID
6432    Synechococcus sp. JA-3-3Ab
80    Synechococcus sp. JA-2-3B'a(2-13)
19    synthetic construct

tabchk

A huge chunk of my time is spent doing ETL operations -- extract, transform, load -- meaning someone sends me data (Excel or delimited-text, JSON/XML), and I put it into some sort of database. I usually want to inspect the data to see what it looks like, and it's hard to see the data when it's in columnar format. I'd rather see it formatted vertically. E.g., here is an export of some iMicrobe sample annotations:

$ tabchk.py -d imicrobe-blast-annots.txt
// ****** Record 1 ****** //
sample_id             : 1
project_name          : Acid Mine Drainage Metagenome
ontologies            : ENVO:01000033 (oceanic pelagic zone biome),ENVO:00002149 (sea water),ENVO:00000015 (ocean)
project_id            : 1
source_mat_id         : 33
sample_type           : metagenome
sample_description    : ACID_MINE_02 - 5-Way (CG) Acid Mine Drainage Biofilm
sample_name           : ACID_MINE_02
collection_start_time : 2002-03-01 00:00:00.0
collection_stop_time  : 2002-03-01 00:00:00.0
site_name             : Richmond Mine
latitude              : 40.666668
longitude             : -122.51667
site_description      : Acid Mine Drainage
country               : UNITED STATES
region                : Iron Mountain, CA
library_acc           : JGI_AMD_5WAY_IRNMTN_LIB_20020301
sequencing_method     : dideoxysequencing (Sanger)
dna_type              : gDNA
num_of_reads          : 180713
habitat_name          : waste water
temperature           : 42
sample_acc            : JGI_AMD_5WAY_IRNMTN_SMPL_20020301

If you open that file in a text editor, you will see there are many more column headers than are shown here. The lines are extremely long, and it's a pretty sparse matrix of data. The -d flag to the program indicates to show a "dense" matrix, i.e., leave out the empty fields. I find the above format much easier to read. Here is the program to do that:

$ cat -n tabchk.py
     1    #!/usr/bin/env python3
     2
     3    import argparse
     4    import csv
     5    import os
     6    import sys
     7
     8    # --------------------------------------------------
     9    def get_args():
    10        parser = argparse.ArgumentParser(description='Check a delimited text file')
    11        parser.add_argument('file', metavar='str', help='File')
    12        parser.add_argument('-s', '--sep', help='Field separator',
    13                            metavar='str', type=str, default='\t')
    14        parser.add_argument('-l', '--limit', help='How many records to show',
    15                            metavar='int', type=int, default=1)
    16        parser.add_argument('-d', '--dense', help='Not sparse (skip empty fields)',
    17                            action='store_true')
    18        return parser.parse_args()
    19
    20    # --------------------------------------------------
    21    def main():
    22        args = get_args()
    23        file = args.file
    24        limit = args.limit
    25        sep = args.sep
    26        dense = args.dense
    27
    28        if not os.path.isfile(file):
    29            print('"{}" is not a file'.format(file))
    30            sys.exit(1)
    31
    32        with open(file) as csvfile:
    33            reader = csv.DictReader(csvfile, delimiter=sep)
    34
    35            for i, row in enumerate(reader):
    36                vals = dict([x for x in row.items() if x[1] != '']) if dense else row
    37                flds = vals.keys()
    38                longest = max(map(len, flds))
    39                fmt = '{:' + str(longest + 1) + '}: {}'
    40                print('// ****** Record {} ****** //'.format(i+1))
    41                for key, val in vals.items():
    42                    print(fmt.format(key, val))
    43
    44                if i + 1 == limit:
    45                    break
    46
    47    # --------------------------------------------------
    48    if __name__ == '__main__':
    49        main()

FASTA

Now let's finally get into parsing good, old FASTA files. The argparse and csv modules are standard in Python, but for FASTA we're going to need to install the BioPython (http://biopython.org/\ module. This should work for you:

$ python3 -m pip install biopython

For this exercise, I'll use a few reads from the Global Ocean Sampling Expedition (https://imicrobe.us/sample/view/578\:

$ head -5 CAM_SMPL_GS108.fa
>CAM_READ_0231669761 /library_id="CAM_LIB_GOS108XLRVAL-4F-1-400" /sample_id="CAM_SMPL_GS108" raw_id=SRA_ID=SRR066139.70645 raw_id=FG67BMZ02PUFIF
ATTTACAATAATTTAATAAAATTAACTAGAAATAAAATATTGTATGAAAATATGTTAAATAATGAAAGTTTTT
CAGATCGTTTAATAATATTTTTCTTCCATTTTGCTTTTTTCTAAAATTGTTCAAAAACAAACTTCAAAGGAAA
ATCTTCAAAATTTACATGATTTTATATTTAAACAAATAGAGTTAAGTATAAGAGAAATTGGATATGGTGATGC
TTCAATAAATAAAAAAATGAAAGAGTATGTCAATGTGATGTACGCAATAATTGACAAAGTTGATTCATGGGAA

Let's write a script that mimics seqmagick:

$ seqmagick info CAM_SMPL_GS108.fa
name              alignment    min_len   max_len   avg_len  num_seqs
CAM_SMPL_GS108.fa FALSE            168       440    325.60         5

We'll skip the "alignment" and just do min/max/avg lengths and the number of sequences. You can pretty much copy and paste the example code from http://biopython.org/wiki/SeqIO. Here is the output from our script, "seqmagique.py":

$ ./seqmagique.py *.fa
name                 min_len    max_len    avg_len   num_seqs
CAM_SMPL_GS108.fa        168        440      325.6          5
CAM_SMPL_GS112.fa         71        517      398.4          5

The code to produce this builds on our earlier skills of lists and dictionaries as we will parse each file and save a dictionary of stats into a list, then we will iterate over that list at the end to show the output.

$ cat -n seqmagique.py
     1    #!/usr/bin/env python3
     2    """Mimic seqmagick, print stats on FASTA sequences"""
     3
     4    import os
     5    import sys
     6    from statistics import mean
     7    from Bio import SeqIO
     8
     9    files = sys.argv[1:]
    10
    11    if len(files) < 1:
    12        print('Usage: {} F1.fa [F2.fa...]'.format(os.path.basename(sys.argv[0])))
    13        sys.exit(1)
    14
    15    info = []
    16    for file in files:
    17        lengths = []
    18        for record in SeqIO.parse(file, "fasta"):
    19            lengths.append(len(record.seq))
    20
    21        info.append({'name': os.path.basename(file),
    22                     'min_len': min(lengths),
    23                     'max_len': max(lengths),
    24                     'avg_len': mean(lengths),
    25                     'num_seqs': len(lengths)})
    26
    27    if len(info):
    28        longest_file_name = max([len(f['name']) for f in info])
    29        fmt = '{:' + str(longest_file_name) + '} {:>10} {:>10} {:>10} {:>10}'
    30        flds = ['name', 'min_len', 'max_len', 'avg_len', 'num_seqs']
    31        print(fmt.format(*flds))
    32        for rec in info:
    33            print(fmt.format(*[rec[fld] for fld in flds]))
    34    else:
    35        print('I had trouble parsing your data')

FASTA subset

Sometimes you may only want to use part of a FASTA file, e.g., you want the first 1000 sequences to test some code, or you have samples that vary wildly in size and you want to sub-sample them down to an equal number of reads. Here is a Python program that will write the first N samples to a given output directory:

$ cat -n fa_subset.py
     1    #!/usr/bin/env python3
     2    """Subset FASTA files"""
     3
     4    import argparse
     5    import os
     6    import sys
     7    from Bio import SeqIO
     8
     9    # --------------------------------------------------
    10    def get_args():
    11        """get args"""
    12        parser = argparse.ArgumentParser(description='Split FASTA files')
    13        parser.add_argument('fasta', help='FASTA input file', metavar='FILE')
    14        parser.add_argument('-n', '--num', help='Number of records per file',
    15                            type=int, metavar='NUM', default=500000)
    16        parser.add_argument('-o', '--out_dir', help='Output directory',
    17                            type=str, metavar='DIR', default='subset')
    18        return parser.parse_args()
    19
    20    # --------------------------------------------------
    21    def main():
    22        """main"""
    23        args = get_args()
    24        fasta = args.fasta
    25        out_dir = args.out_dir
    26        num_seqs = args.num
    27
    28        if not os.path.isfile(fasta):
    29            print('--fasta "{}" is not valid'.format(fasta))
    30            sys.exit(1)
    31
    32        if os.path.dirname(fasta) == out_dir:
    33            print('--outdir cannot be the same as input files')
    34            sys.exit(1)
    35
    36        if num_seqs < 1:
    37            print("--num cannot be less than one")
    38            sys.exit(1)
    39
    40        if not os.path.isdir(out_dir):
    41            os.mkdir(out_dir)
    42
    43        basename = os.path.basename(fasta)
    44        out_file = os.path.join(out_dir, basename)
    45        out_fh = open(out_file, 'wt')
    46        num_written = 0
    47
    48        for record in SeqIO.parse(fasta, "fasta"):
    49            SeqIO.write(record, out_fh, "fasta")
    50            num_written += 1
    51
    52            if num_written == num_seqs:
    53                break
    54
    55        print('Done, wrote {} sequence{} to "{}"'.format(
    56            num_written, '' if num_written == 1 else 's', out_file))
    57
    58    # --------------------------------------------------
    59    if __name__ == '__main__':
    60        main()

FASTA splitter

I seem to have implemented my own FASTA splitter a few times in as many languages. Here is one that writes a maximum number of sequences to each output file. It would not be hard to instead write a maximum number of bytes, but, for the short reads I usually handle, this works fine. Again I will use the BioPython SeqIO module to parse the FASTA files

$ cat -n fasplit.py
     1    #!/usr/bin/env python3
     2    """split FASTA files"""
     3
     4    import argparse
     5    import os
     6    from Bio import SeqIO
     7
     8    # --------------------------------------------------
     9    def main():
    10        """main"""
    11        args = get_args()
    12        fasta = args.fasta
    13        out_dir = args.out_dir
    14        max_per = args.num
    15
    16        if not os.path.isfile(fasta):
    17            print('--fasta "{}" is not valid'.format(fasta))
    18            exit(1)
    19
    20        if not os.path.isdir(out_dir):
    21            os.mkdir(out_dir)
    22
    23        if max_per < 1:
    24            print("--num cannot be less than one")
    25            exit(1)
    26
    27        i = 0
    28        nseq = 0
    29        nfile = 0
    30        out_fh = None
    31        basename, ext = os.path.splitext(os.path.basename(fasta))
    32
    33        for record in SeqIO.parse(fasta, "fasta"):
    34            if i == max_per:
    35                i = 0
    36                if out_fh is not None:
    37                    out_fh.close()
    38                    out_fh = None
    39
    40            i += 1
    41            nseq += 1
    42            if out_fh is None:
    43                nfile += 1
    44                path = os.path.join(out_dir, basename + '.' + str(nfile) + ext)
    45                out_fh = open(path, 'wt')
    46
    47            SeqIO.write(record, out_fh, "fasta")
    48
    49        print('Done, wrote {} sequence{} to {} file{}'.format(
    50            nseq, '' if nseq == 1 else 's',
    51            nfile, '' if nfile == 1 else 's'))
    52
    53    # --------------------------------------------------
    54    def get_args():
    55        """get args"""
    56        parser = argparse.ArgumentParser(description='Split FASTA files')
    57        parser.add_argument('-f', '--fasta', help='FASTA input file',
    58                            type=str, metavar='FILE', required=True)
    59        parser.add_argument('-n', '--num', help='Number of records per file',
    60                            type=int, metavar='NUM', default=50)
    61        parser.add_argument('-o', '--out_dir', help='Output directory',
    62                            type=str, metavar='DIR', default='fasplit')
    63        return parser.parse_args()
    64
    65    # --------------------------------------------------
    66    if __name__ == '__main__':
    67        main()

If you type make in the "python/fasta-splitter" directory, you should see:

$ make
./fasplit.py -f POV_L.Sum.O.1000m_reads.fa -o pov -n 100
Done, wrote 2061 sequences to 21 files

We can verify that things worked:

$ for file in *; do echo -n $file && grep '^>' $file | wc -l; done
POV_L.Sum.O.1000m_reads.1.fa     100
POV_L.Sum.O.1000m_reads.10.fa     100
POV_L.Sum.O.1000m_reads.11.fa     100
POV_L.Sum.O.1000m_reads.12.fa     100
POV_L.Sum.O.1000m_reads.13.fa     100
POV_L.Sum.O.1000m_reads.14.fa     100
POV_L.Sum.O.1000m_reads.15.fa     100
POV_L.Sum.O.1000m_reads.16.fa     100
POV_L.Sum.O.1000m_reads.17.fa     100
POV_L.Sum.O.1000m_reads.18.fa     100
POV_L.Sum.O.1000m_reads.19.fa     100
POV_L.Sum.O.1000m_reads.2.fa     100
POV_L.Sum.O.1000m_reads.20.fa     100
POV_L.Sum.O.1000m_reads.21.fa      61
POV_L.Sum.O.1000m_reads.3.fa     100
POV_L.Sum.O.1000m_reads.4.fa     100
POV_L.Sum.O.1000m_reads.5.fa     100
POV_L.Sum.O.1000m_reads.6.fa     100
POV_L.Sum.O.1000m_reads.7.fa     100
POV_L.Sum.O.1000m_reads.8.fa     100
POV_L.Sum.O.1000m_reads.9.fa     100

GFF

Two of the most common output files in bioinformatics, GFF (General Feature Format) and BLAST's tab/CSV files do not include headers, so it's up to you to merge in the headers. Additionally, some of the lines may be comments (they start with "#" just like bash and Python), so you should skip those. Further, the last field in GFF is basically a dumping ground for whatever else the data provider felt like putting there. Usually it's a bunch of "key=value" pairs, but there's no guarantee. Let's take a look at parsing the GFF output from Prodigal:

$ cat -n parse_prodigal_gff.py
     1    #!/usr/bin/env python3
     2
     3    import argparse
     4    import csv
     5    import os
     6    import sys
     7
     8    # --------------------------------------------------
     9    def get_args():
    10        parser = argparse.ArgumentParser(description='Parse Prodigal GFF')
    11        parser.add_argument('gff', metavar='GFF', help='Prodigal GFF output')
    12        parser.add_argument('-m', '--min', help='Minimum score',
    13                            metavar='float', type=float, default=0.0)
    14        return parser.parse_args()
    15
    16    # --------------------------------------------------
    17    def main():
    18        args = get_args()
    19        gff_file = args.gff
    20        min_score = args.min
    21        flds = 'seqname source feature start end score strand frame attribute'.split()
    22
    23        for line in open(gff_file):
    24            if line[0] == '#':
    25                continue
    26
    27            row = dict(zip(flds, line.split('\t')))
    28            attrs = {}
    29            if 'attribute' in row:
    30                for fld in row['attribute'].split(';'):
    31                    if '=' in fld:
    32                        name, val = fld.split('=')
    33                        attrs[name] = val
    34
    35            if 'score' in attrs and float(attrs['score']) > min_score:
    36                print(row)
    37                break
    38
    39    # --------------------------------------------------
    40    if __name__ == '__main__':
    41        main()

results matching ""

    No results matching ""