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

写在前面

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

image.png

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

使用Argparse模块

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可能比较合适。