写在前面
在上一则推文《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
得到这个提示
image.png
运行
python testArgsParse.py --inSam Treat.20M.merged.sam --outSam dedup.sorted.sam --maxHitPos 3
没问题，maxHitPos一般设置为20可能比较合适。
网友评论