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

写在前面

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

《用python做生物信息数据分析(3-处理命令行参数)》 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

得到这个提示

《用python做生物信息数据分析(3-处理命令行参数)》 image.png

运行

python testArgsParse.py --inSam Treat.20M.merged.sam --outSam dedup.sorted.sam --maxHitPos 3

没问题,maxHitPos一般设置为20可能比较合适。

    原文作者:生信札记
    原文地址: https://www.jianshu.com/p/148e28f141db
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞