用python做生物信息数据分析（3-处理命令行参数）

写在前面

在上一则推文《2-pysam》中，写了一个基本能用的脚本。但是，很明显，那个可用于生产的脚本。因为他不能很好的处理命令行参数。总的来说，就是一个死的脚本，见下图

image.png

其中输入文件Treat.20M.sorted.sam写死在脚本中。如果下一次有新的文件要处理，那么就必须修改脚本。要让他脚本活过来，就需要处理命令行参数。在perl里面，我直接是从零码起，再做判断。在java上，我自己写了一个ArgsParser的类，自以为还不错。而在python里面，经过基本检索，可以确定，应该是使用Argparse

使用Argparse模块

查看一下模块的文档，
https://docs.python.org/3/howto/argparse.html
目标很简单，就是能接受三个参数
稍微看了下文档，可能我们需要的是

import argparse
parser = argparse.ArgumentParser()
parser.add_argument("square", help="display a square of a given number",
                    type=int)
args = parser.parse_args()

其中的type应该是可选项，随后就可以做自己想做的事情

脚本调整

原来的脚本

import pysam
# filter sam file to remove reads mapped to repeat regions.\
samfile = pysam.AlignmentFile("Treat.20M.merged.sam")
tmpfile = pysam.AlignmentFile("dedup.sam", "w", template=samfile)
lineCount = 0
max_hit_num = 3
pre_read_id = ""
cur_read_list = list()
for read in samfile:
    cur_read_id = read.qname
    if cur_read_id == pre_read_id:
        cur_read_list.append(read)
    else:
        if len(cur_read_list) < max_hit_num:
            for cur_read in cur_read_list:
                tmpfile.write(cur_read)
        cur_read_list.clear()
        cur_read_list.append(read)
        pre_read_id = cur_read_id
    lineCount = lineCount + 1
    if lineCount > 100:
        break
if len(cur_read_list) != 0 & len(cur_read_list) < max_hit_num:
    for cur_read in cur_read_list:
        tmpfile.write(cur_read)
cur_read_list.clear()
# sort sam file
pysam.sort("-o", "dedup.sorted.sam", "dedup.sam")

抽取需要设置的参数，整体获得三个

import pysam
in_sam = "Treat.20M.merged.sam"
out_sam = "dedup.sorted.sam"
max_hit_num = 3
# 
tmp_file = in_sam + ".dedup.sam";
# filter sam file to remove reads mapped to repeat regions.
samfile = pysam.AlignmentFile(in_sam)
tmpfile = pysam.AlignmentFile(tmp_file, "w", template=samfile)
# lineCount = 0
pre_read_id = ""
cur_read_list = list()
for read in samfile:
    cur_read_id = read.qname
    if cur_read_id == pre_read_id:
        cur_read_list.append(read)
    else:
        if len(cur_read_list) < max_hit_num:
            for cur_read in cur_read_list:
                tmpfile.write(cur_read)
        cur_read_list.clear()
        cur_read_list.append(read)
        pre_read_id = cur_read_id
    # lineCount = lineCount + 1
    # if lineCount > 100:
    #    break
if len(cur_read_list) != 0 & len(cur_read_list) < max_hit_num:
    for cur_read in cur_read_list:
        tmpfile.write(cur_read)
cur_read_list.clear()
# sort sam file
pysam.sort("-o",out_sam , tmp_file)

使用ArgsParse调整脚本

import pysam
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--inSam", help="set input sam file, should be name-sorted, single-end reads")
parser.add_argument("--outSam", help="set output sam file")
parser.add_argument("--maxHitPos", help="set max hit pos to define repeat-generated read",type=int)
args = parser.parse_args()

in_sam = args.inSam
out_sam = args.outSam
max_hit_num = args.maxHitPos
# in_sam = "Treat.20M.merged.sam"
# out_sam = "dedup.sorted.sam"
# max_hit_num = 3
#
tmp_file = in_sam + ".dedup.sam";
# filter sam file to remove reads mapped to repeat regions.
samfile = pysam.AlignmentFile(in_sam)
tmpfile = pysam.AlignmentFile(tmp_file, "w", template=samfile)
# lineCount = 0
pre_read_id = ""
cur_read_list = list()
for read in samfile:
    cur_read_id = read.qname
    if cur_read_id == pre_read_id:
        cur_read_list.append(read)
    else:
        if len(cur_read_list) < max_hit_num:
            for cur_read in cur_read_list:
                tmpfile.write(cur_read)
        cur_read_list.clear()
        cur_read_list.append(read)
        pre_read_id = cur_read_id
    # lineCount = lineCount + 1
    # if lineCount > 100:
    #    break
if len(cur_read_list) != 0 & len(cur_read_list) < max_hit_num:
    for cur_read in cur_read_list:
        tmpfile.write(cur_read)
cur_read_list.clear()
samfile.close()
tmpfile.close()
# sort sam file
pysam.sort("-o",out_sam , tmp_file)

使用脚本

查看help。注意，如果参数名对应的字符串是inSam，那么就是一个位置参数，这个其实不是太好用，而加上--，就变成--inSam，则可以不错是有位置随意的参数，所以我都加上了

python testArgsParse.py -h

得到这个提示