12下一页
返回列表 发新帖

【python练习题二】正则表达式

[复制链接]

5

主题

20

帖子

318

积分

版主

Rank: 7Rank: 7Rank: 7

积分
318
发表于 2019-11-15 11:28:28 | 显示全部楼层 | 阅读模式
本帖最后由 chenming 于 2019-11-15 11:33 编辑

正则表达式即为字符串的模式匹配,举个例子,在草原上的马群中挑选出长角的白色马匹,这里的长角/白色/马都是限定词,也就是我们用来构建一个“模式”的条件,
正则表达式中有很多抽象的符号,这些符号都是为了满足我们各种抽象而特别的要求而设计的,为了记住这些符号,我们可以先对这些特殊符号进行了解,再慢慢通过例子
和练习加深印象。
python处理正则表达式的模块叫做re模块,请试着回答以下问题:(不要求每个人都能回答所有问题,只需要根据自己所知道的,在网络上所查到的,感兴趣的内容进行回答,相互补充,每个人都可以获得正则表达式的全貌,习题请尽量完成。)
  • 说出re模块最常用的7个函数或方法,并说明它们的区别和各自的应用场景
  • 请试着说出以下的特殊字符的含义
  1. ^
  2. $
  3. .
  4. *
  5. [a-z]
  6. [^0-9]
  7. {2,4}
  8. +
  9. ?
复制代码
3.请试着说出以下特殊字符的含义
  1. \w \W \B \b \d \D \s \S
复制代码
4.请说明什么是正则表达式中的贪婪与非贪婪匹配
  5.请说明如何让正则表达式可以忽略大小写字符/进行多行匹配?
  6.完成以下习题:
(1)使用re模块从形如以下截图的vcf文件中获取INFO列的基因名称
文件请使用/WORK/Disease/chenming/Workshop/re_exercise/re_exercise.vcf.gz

(2)请使用re模块找到人类GRCh37参考基因组22号染色体中所有AATT模式的起始和终止位置
(3)请使用re模块找出两个样本的共有和各自独有的变异位点,这里认为./.和0/0均为未突变,请
使用文件/WORK/Disease/chenming/Workshop/re_exercise/snp.indel.merged.freq.func.syn.deleterious.xls

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x

0

主题

1

帖子

36

积分

新手上路

Rank: 1

积分
36
发表于 2019-11-15 18:44:14 | 显示全部楼层
本帖最后由 张子龙 于 2019-11-15 18:46 编辑

1、说出re模块最常用的7个函数或方法,并说明它们的区别和各自的应用场景
re.search 在整个字符串中查找,返回第一个匹配内容,如果找到返回match对象,没找到返回None
re.match 只匹配字符串的开始,返回第一个匹配内容,如果字符串不符合正则表达式,则匹配失败,函数返回None
re.split 按照匹配的字符串进行分割
re.sub 替换匹配的字符串,返回替换后的字符串
re.findall 在字符串中找到正则表达式所匹配的所有字符串,并返回一个列表,如果没有找到匹配的,则返回空列表
re.compile 编译正则表达式,生成一个正则表达式对象,供match和search这两个函数使用
re.finditer 在字符串中找到正则表达式所匹配的所有子串,并把它们作为一个迭代器返回
2、请试着说出以下的特殊字符的含义
^    从第一个字符串开始匹配
$    字符串的末尾
.    除换行符以外的任意字符
*    数目>=0
[a-z]    匹配小写字母
[^0-9]    匹配非数字
{2,4}    模式出现2-4次
+    数目>=1
?    出现或者不出现
3.请试着说出以下特殊字符的含义
\w 字母数字下划线
\W 非字母数字下划线
\B 非单词边界
\b 单词边界
\d 数字
\D 非数字
\s 空白符
\S 非空白符
4.请说明什么是正则表达式中的贪婪与非贪婪匹配
贪婪匹配:尽量匹配更多的满足条件的字符
非贪婪匹配:尽量匹配最少的满足条件的字符
5.请说明如何让正则表达式可以忽略大小写字符/进行多行匹配?
忽略大小写:re.compile(,re.I)
多行匹配:re.compile(,re.M)
6.完成以下习题:
(1)
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. __metaclass__ = type

  4. import re
  5. import gzip
  6. fr_file = gzip.open("/WORK/Disease/chenming/Workshop/re_exercise/re_exercise.vcf.gz")
  7. for line in fr_file:
  8.     re_result = re.search(r"GENE=(.+?);", line)
  9.     if re_result:
  10.         print re_result.group(1)
复制代码
(2)
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. __metaclass__ = type

  4. import re
  5. human_ref = file("/PUBLIC/database/HUMAN/genome/Human/human_g1k_v37_decoy.fasta")
  6. str22_sequence = ""
  7. for line in human_ref:
  8.     if line.startswith(">"):
  9.         if line.startswith(">22"):
  10.             flag = True
  11.             continue
  12.         else:
  13.             flag = False
  14.             continue
  15.     if flag:
  16.         str22_sequence += line.strip()

  17. for match in re.finditer(r"AATT",str22_sequence):
  18.     print match.span(),str22_sequence[match.start():match.end()]
  19.         
复制代码
(3)
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3. __metaclass__ = type

  4. import re
  5. re_not_have_mutation = re.compile(r"(\./\.)|(0/0)")
  6. fr_file = file("snp.indel.merged.freq.func.syn.deleterious.xls")
  7. fr_file.next()
  8. for line in fr_file:
  9.     line_list = line.strip("\n").split("\t")
  10.     s1_format = line_list[55]
  11.     s2_format = line_list[56]
  12.     if not re.search(re_not_have_mutation, s1_format) and re.search(re_not_have_mutation, s2_format):
  13.         print "S1单独有的突变",s1_format,s2_format
  14.     elif  re.search(re_not_have_mutation, s1_format) and not re.search(re_not_have_mutation, s2_format):
  15.         print "S2单独有的突变",s1_format,s2_format
  16.     elif not re.search(re_not_have_mutation, s1_format) and not re.search(re_not_have_mutation, s2_format):
  17.         print "共有的突变",s1_format,s2_format
复制代码



回复

使用道具 举报

2

主题

6

帖子

171

积分

注册会员

Rank: 2

积分
171
发表于 2019-11-18 16:33:10 | 显示全部楼层
本帖最后由 lmt 于 2019-11-18 16:50 编辑

问题一:

"""re 函数
"""
# compile(pattern, flags=0)           # 使用任何可选的标记来编译正则表达式的模式,然后返回一个正则表达式对象
# match(pattern, string, flags=0)     # 尝试使用带有可选的标记的正则表达式的模式来匹配字符串。如果匹配成功,就返回匹配对象; 如果失败,就返回 None
# search(pattern, string, flags=0)    # 使用可选标记搜索字符串中第一次出现的正则表达式模式。 如果匹配成功,则返回匹配对象; 如果失败,则返回 None
# findall(pattern, string, flags=0)   # 查找字符串中所有(非重复)出现的正则表达式模式,并返回一个匹配列表
# finditer(pattern, string, flags=0)  # 与 findall()函数相同,但返回的不是一个列表,而是一个迭代器。 对于每一次匹配,迭代器都返回一个匹配对象
# split(pattern, string, max=0)       # 根据正则表达式的模式分隔符, split函数将字符串分割为列表,然后返回成功匹配的列表,分隔最多操作 max 次(默认分割所有匹配成功的位置)
# sub(pattern, repl, string, count)   # 使用 repl 替换所有正则表达式的模式在字符串中出现的位置,除非定义 count, 否则就将替换所有出现的位置( 另见 subn()函数,该函数返回替换操作的数目)
# purge()                             # 清除隐式编译的正则表达式模式



问题二:

^:匹配行首
$:匹配行尾
.:匹配任意字符(除了\n)
*:匹配前一个字符出现0次或者无限次,即可有可无
[a-z]: 中括号表示匹配[]中列举的字符;实例表示匹配a到z的小写字母
[^0-9]:[^]方括号总的^表示取反的意思,实例表示不匹配0到9的数字
{2,4}:匹配前一个字符出现2至4次
+:匹配前一个字符出现1次或者无限次,至少出现1次
?:匹配前面一个字符出现1次或者0次,即最多1次


问题三:

  1. \w:匹配单词字符
  2. \W:匹配非单词字符

  3. re.match('\w', '_a')
  4. #  <re.Match object; span=(0, 1), match='_'>
复制代码
  1. \b:表示字母数字与非字母数字的边界,     非字母数字与字母数字的边界。

  2. \B:表示字母数字与(非非)字母数字的边界,非字母数字与非字母数字的边界。

  3. re.split('123\\b', '==123!! abc123. 123. 123abc. 123')
  4. # ['==', '!! abc', '. ', '. 123abc. ', '']
复制代码
  1. \d:匹配数字。 即0-9
  2. \D:匹配非数字,即不是数字

  3. re.match('\d\d', '123')
  4. # <re.Match object; span=(0, 2), match='12'>
复制代码
  1. \s:匹配空格, tab键, \n
  2. \S:匹配非空格

  3. re.match('\s', ' a')
  4. # <re.Match object; span=(0, 1), match=' '>

  5. re.match('\S',  ' a')
  6. # None

复制代码
问题四:贪婪匹配 和 非贪婪匹配
D:/youdaonote/m13554221497@163.com/9c7e81729c12420fa69348c367b3afa2/clipboard.png

贪婪匹配:尽可能匹配多的字符串
非贪婪匹配:匹配到结果就好,总是尝试尽可能少的字符。(比较佛系)

  1. import re
  2. match_str = 'abcbcedfcdggc'
  3. pattern1 = r'ab.*c'
  4. pattern2 = r'ab.*?c'

  5. # 贪婪匹配
  6. re.match(pattern1, match_str)
  7. # <re.Match object; span=(0, 13), match='abcbcedfcdggc'>

  8. # 非贪婪匹配
  9. re.match(pattern2, match_str)
  10. # <re.Match object; span=(0, 3), match='abc'>
复制代码
问题五:正则的compile函数
思想:先将需要匹配的字符串,写成pattern实例;然后使用pattern实例处理文本并获取匹配结果。
用法: re.compile(pattern, flags=0)

  1. pattern : 一个字符串形式的正则表达式
  2. flags 可选,表示匹配模式,比如忽略大小写,多行模式等,具体参数为:
  3. <b>re.I 忽略大小写</b>
  4. re.L 表示特殊字符集 \w, \W, \b, \B, \s, \S 依赖于当前环境
  5. <b>re.M 多行模式</b>
  6. re.S 即为' . '并且包括换行符在内的任意字符(' . '不包括换行符)
  7. re.U 表示特殊字符集 \w, \W, \b, \B, \d, \D, \s, \S 依赖于 Unicode 字符属性数据库
  8. re.X 为了增加可读性,忽略空格和' # '后面的注释
复制代码
练习:

练习一:

  1. #!/usr/bin/env python
  2. #-*- coding:utf-8 -*-
  3. # get gene name
  4. import re
  5. import sys
  6. def get_genename(infile):
  7.     gene_list = []
  8.     if infile.endswith('gz'):
  9.         import gzip
  10.         f = gzip.open(infile, 'r')
  11.     else:
  12.         f = open(infile, 'r')
  13.    
  14.     for line in f:
  15.         if line.startswith('#'):
  16.             pass
  17.         else:
  18.             result = re.search(r'GENE=(.*?);', line)
  19.             if result:
  20.                 genename = result.group(1)
  21.                 gene_list.append(genename)

  22.     return list(set(gene_list))

  23. if __name__ == "__main__":
  24.     infile = sys.argv[1]
  25.     get_genename(infile)
复制代码
  • 练习三
  1. #!/usr/bin/env python
  2. #-*- coding:utf-8 -*-
  3. import re
  4. import sys
  5. def check_genetype(infile):
  6.     with open(infile, 'r')  as indata:
  7.         header = indata.readline().strip().split('\t')
  8.         s1_pos = header.index('S1')
  9.         s2_pos = header.index('S2')
  10.         pattern_no_mut = re.compile(r'(\./\.)|(0/0)')
  11.         for line in indata:
  12.             lline = line.strip().split('\t')
  13.             s1_info = lline[s1_pos]
  14.             s2_info = lline[s2_pos]
  15.            
  16.             s1_tp = re.search(pattern_no_mut, s1_info)
  17.             s2_tp = re.search(pattern_no_mut, s2_info)
  18.             if not s1_tp and not s2_tp:
  19.                 print('共有变异位点: %s\t%s' % (s1_info, s2_info))
  20.             elif not s1_tp and s2_tp:
  21.                 print('just s1 have : %s\t%s' %  (s1_info, s2_info))
  22.             elif s1_tp and not s2_tp:
  23.                 print('just s2 have: %s\t%s ' %(s1_info, s2_info) )
  24.             else:
  25.                 print('都没有的突变 %s\t%s' % (s1_info, s2_info))

  26.    
  27. if __name__ == "__main__":
  28.     infile = sys.argv[1]
  29.     check_genetype(infile)
复制代码



   

D:/youdaonote/m13554221497@163.com/9c7e81729c12420fa69348c367b3afa2/clipboard.png



回复

使用道具 举报

5

主题

20

帖子

318

积分

版主

Rank: 7Rank: 7Rank: 7

积分
318
发表于 2019-11-19 10:23:46 | 显示全部楼层
子龙和梦婷对于re模块的主要内容介绍的非常详细,很棒!
回复

使用道具 举报

5

主题

20

帖子

318

积分

版主

Rank: 7Rank: 7Rank: 7

积分
318
发表于 2019-11-19 10:37:00 | 显示全部楼层
练习题的话,最好能对得到的结果进行一个简单的描写,例如练习题一,你找到了多少个基因,练习题二,你找到了多少个AATT,练习题三,你找到了多少共有和各自特有的位点
回复

使用道具 举报

6

主题

11

帖子

576

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
576

优秀版主荣誉管理论坛元老

发表于 2019-11-20 14:54:35 | 显示全部楼层
习题3
  • 输入文件可以是文本,也可以是gzip压缩文件
  • 可以统计2个样本或多个样本的特有,共有
  • 输出结果文件,打印统计信息

源码:
https://gitee.com/suqingdong/codes/ryja94tqdsfu7ko1nehxg15

运行示例:



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

6

主题

11

帖子

576

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
576

优秀版主荣誉管理论坛元老

发表于 2019-11-20 15:44:01 | 显示全部楼层
习题2
  • 可指定参考基因组,染色体,查询序列
  • 使用pysam快速定位染色体
  • 每次读取部分数据,减少内存占用
  • 重复匹配次数较多时,使用re.compile预编译,以提升性能


源码:
https://gitee.com/suqingdong/codes/p9siwd52y8bva7fzgutjn17

运行示例:

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

0

主题

1

帖子

6

积分

新手上路

Rank: 1

积分
6
发表于 2019-12-5 23:03:37 | 显示全部楼层
1、re模块最常用的7个函数
    match():从头匹配,从字符串的第一个字符开始匹配,如果找到返回匹配的对象,没找到返回None
    search():从整个文本搜索,在整个字符串中查找,返回第一个匹配内容,如果找到返回查找的对象,没找到返回None
    findall():找到所有符合的,在字符串中匹配,如果成功返回match对象,如果失败返回None
    split():分割,按照匹配的字符串进行分割
    sub():替换,替换匹配的子字符串,返回替换之后的字符串


2、特殊字符的含义
    ^:匹配字符串开头
    $:匹配字符串结尾
    . :匹配换行符外的任意字符
    *:匹配前面出现的正则表达式零次或多次
    [a-z]:匹配从a到z的任意一个字符
    [^0-9]:不匹配0-9中的任意一个字符
    {2,4}:匹配前面重复2-4次的正则表达式
    +:匹配前面出现的正则表达式一次或多次
?:匹配前面出现的正则表达式零次或一次
\w:匹配任意数字字母字符
\W:匹配任意非数字字母字符
\B:匹配非字母和数字边界
\b:匹配字母或数字边界
\d:匹配任意数字
\D:匹配任意非数字字符
\s:匹配任意空白符
\S:匹配任意非空白符

3、贪婪匹配趋向于匹配最大长度,会返回最长的匹配对象;非贪婪匹配趋向于匹配尽量少的结果,一旦符合匹配条件,就停止匹配。

4、忽略大小写:flags=re.I
   多行匹配:flags=re.M



练习题一
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import re
import argparse
import os

parser=argparse.ArgumentParser(description="extract site gene info")
parser.add_argument('--inf',help="The path of input file",required=True)
parser.add_argument('--out',help="The path of output file",required=True)

args = parser.parse_args()

infile= args.inf
out= args.out

def safe_open(file,mode):
    file = os.path.abspath(file)
    if not os.path.exists(file):
        exit('%s is not exists.' %file)
    elif file.endswith('.gz'):
        import gzip
        return gzip.open(file,mode)
    else:
        return open(file,mode)

with safe_open(infile,'r') as inf,open(out,'w') as outf:
    for line in inf:
       if line.startswith('#'):
           pass
       else:
           target_word=re.search(r'(.*)GENE=(.*?);.*',line)
          # print target_word.group(2)
           outf.write(target_word.group(2)+'\n')
回复

使用道具 举报

1

主题

2

帖子

23

积分

新手上路

Rank: 1

积分
23
发表于 2019-12-6 11:50:54 | 显示全部楼层
1.re模块最常用的7个函数和方法,各自的区别和应用场景。
  (1)re.match 尝试从字符串的起始位置匹配一个模式,如果不是起始位置匹配成功的话,match就返回none。
    用法:re.match(pattern,string,flags=0)
    参数说明:pattern 匹配的正则表达式;string 要匹配的字符串;flags 标志位,用于控制正则表达式的匹配方式,如:是否区分大小写,多行匹配等等。
  (2)re.search 扫描整个字符串并返回第一个成功的匹配。
    用法:re.seach(pattern,string,flags=0)
    re.match与re.search的区别:re.match只匹配字符串的开始,如果字符串开始不符合正则表达式,则匹配失败,而re.search匹配整个字符串,直到找到一个匹配。
  (3)re.sub 用于替换字符串中的匹配项。
    用法:re.sub(pattern,repl,string,count=0,flages=0)
    参数:pattern 正则表达式中模式字符串;repl 替换的字符串,也可为一个函数;string 要被查找和替换的原始字符串;count 模式匹配后替换的最大次数,默认0表示替换所有的匹配。
  (4)re.compile 用于编译正则表达式,生成一个正则表达式(pattern)对象,供match()和search()这两个函数使用。
    用法:re.compile(pattern[,flags])
    参数:pattern 一个字符串形式的正则表达式;
         flags 可选,表示匹配模式,比如忽略大小写,多行模式等,具体参数为:re.I忽略大小写,re.L表示特殊字符集\w,\W,\b,\B,\d,\D,\s,\S依赖于当前环境,re.M多行模式,re.S为.并且包括换行符在内的任意字符(.不包括换行符),re.U表示特殊字符集 \w, \W, \b, \B, \d, \D, \s, \S 依赖于 Unicode 字符属性数据库,re.X 为了增加可读性,忽略空格和 # 后面的注释
  (5)findall 在字符串中找到正则表达式所匹配的所有子串,并返回一个列表,如果没有找到匹配的,则返回空列表。
    用法:findall(string[,pos[,endpos]])
    参数:string 待匹配的字符串;pos 可选参数,指定字符串的起始位置,默认为0;endpos 可选参数,指定字符串的结束位置,默认字符串的长度。
    区别:match和search是匹配一次,findall匹配所有。
  (6)re.finditer 和findall类似,在字符串中找到正则表达式中所匹配的所有子串,并把它们作为一个迭代器返回。
    用法:re.finditer(pattern,string,flags=0)
    参数:pattern 匹配正则表达式;string 要匹配的字符串;flags 标志位,用于控制正则表达式的匹配方式,如:是否区分大小写,多行匹配等等。
  (7)re.split 按照能够匹配的子串将字符串分割后返回列表。
    用法:re.split(pattern,string[,maxsplit=0,flags=0])
    参数:pattern 匹配正则表达式;string 要匹配的字符串;flags 标志位,用于控制正则表达式的匹配方式,如:是否区分大小写,多行匹配等等;maxsplit 分隔次数,maxsplit=1分隔1次,默认为0,不限制次数。
2.特殊字符含义
  ^             匹配字符串的开头
  $             匹配字符串的末尾
  .             匹配任意字符,除了换行符,当re.DOTALL标记被指定时,则可以匹配包括换行符的任意字符
  *             匹配0个或多个的表达式
  [a-z]         匹配26个英文小写字母中的任意一个
  [^0-9]        匹配以数字开头的的表达式
  {2,4}         匹配 2 到 4 次由前面的正则表达式定义的片段,贪婪方式
  +             匹配1个或多个的表达式
  ?             匹配0个或1个由前面的正则表达式定义的片段,非贪婪方式
3.\w \W \B \b \d \D \s \S 特殊字符含义
  \w            匹配字母数字及下划线
  \W            匹配非字母数字及下划线
  \B            匹配非单词边界。'er\B' 能匹配 "verb" 中的 'er',但不能匹配 "never" 中的 'er'
  \b            匹配一个单词边界,也就是指单词和空格间的位置。例如, 'er\b' 可以匹配"never" 中的 'er',但不能匹配               "verb" 中的 'er'。
  \d            匹配任意数字,等价于 [0-9].
  \D            匹配任意非数字
  \s            匹配任意空白字符,等价于 [\t\n\r\f]
  \S            匹配任意非空字符
4.正则表达式中的贪婪模式:正则表达式一般趋向于最大长度匹配;非贪婪模式:在整个表达式匹配成功的前提下,尽可能少的匹配。
5.re模块函数的参数flags 表示匹配模式,比如忽略大小写,多行模式等,具体参数为:re.I忽略大小写,re.L表示特殊字符集\w,\W,\b,\B,\d,\D,\s,\S依赖于当前环境,re.M多行模式,re.S为.并且包括换行符在内的任意字符(.不包括换行符),re.U表示特殊字符集 \w, \W, \b, \B, \d, \D, \s, \S 依赖于 Unicode 字符属性数据库,re.X 为了增加可读性,忽略空格和 # 后面的注释。
回复

使用道具 举报

0

主题

4

帖子

128

积分

注册会员

Rank: 2

积分
128
发表于 2019-12-6 12:28:20 | 显示全部楼层
练习1
  1. #!/usr/bin/env python

  2. import re
  3. import sys,gzip

  4. invcf=sys.argv[1]

  5. def openF(infile,mode='r'):
  6.     try:
  7.         if infile.endswith('.gz'):
  8.             mode+='b'
  9.             return gzip.open(infile,mode)
  10.         else:
  11.             return open(infile,mode)
  12.     except Exception as e:
  13.         exit(1)

  14. geneD={}
  15. ptn1=re.compile(r"GENE=(?P<g>.*?);")
  16. #ptn2=re.compile(r'(?<=GENE=)(.+?)(?=;)')
  17. with openF(invcf,'r') as inf:
  18.     for line in inf:
  19.         if line.startswith('#'):
  20.             pass
  21.         else:
  22.             gene=ptn1.search(line).group('g')
  23.             #gene=ptn2.search(line).group(1)
  24.             if gene not in geneD:
  25.                 geneD[gene]=1
  26.             else:
  27.                 geneD[gene]+=1
  28. for k,v in geneD.items():
  29.     print k,":",v
  30. print "Gene Number: ",len(geneD)
复制代码

  1. Gene Number:  115
复制代码


练习2
  1. #!/usr/bin/env python

  2. import re
  3. import sys

  4. infa=sys.argv[1]
  5. pt=sys.argv[2]

  6. ptn1=re.compile(r''+str(pt)+'')
  7. with open(infa,'r') as inf,open(pt+'.xls','w') as outf:
  8.     text=''
  9.     for line in inf:
  10.         if line.startswith('>'):
  11.             pass
  12.         else:
  13.             text+=line.strip()
  14.     n=0
  15.     for seq in re.finditer(ptn1,text):
  16.         n+=1
  17.         outf.write(seq.group()+'\t'+str(seq.span()[0])+'\t'+str(seq.span()[1])+'\n')
  18.     print "find "+pt+" number: ",+n
复制代码

  1. find AATT number:  166912
复制代码


练习3
  1. #!/usr/bin/env python

  2. import re
  3. import sys,gzip

  4. infile=sys.argv[1]

  5. def openF(infile,mode='r'):
  6.     try:
  7.         if infile.endswith('.gz'):
  8.             mode+='b'
  9.             return gzip.open(infile,mode)
  10.         else:
  11.             return open(infile,mode)
  12.     except Exception as e:
  13.         exit(1)

  14. def getTitle(line,t1,t2=''):
  15.     if t1 in line:
  16.         p=line.index(t1)
  17.     elif t2 in line:
  18.         p=line.index(t2)
  19.     else:
  20.         print "There is no ",t1," ",t2,"in the title!"
  21.     return p

  22. def getGT(gt):
  23.     gt1=re.search(r'((?:.|\d))/((?:.|\d))',gt).group(1)
  24.     gt2=re.search(r'((?:.|\d))/((?:.|\d))',gt).group(2)
  25.    
  26.     if gt1==gt2=="0" or gt1==gt2==".":
  27.         flag=0
  28.     elif gt1!=gt2:
  29.         flag=1
  30.     elif gt1 == gt2:
  31.         flag=2
  32.     else:
  33.         print gt
  34.     return flag
  35. share=0
  36. onlys1=0
  37. onlys2=0
  38. diffshare=0
  39. with openF(infile,'r') as inf:
  40.     for line in inf:
  41.         ii=line.strip().split('\t')
  42.         if line.startswith('Priority'):
  43.             s1p=getTitle(ii,'FORMAT')
  44.         else:
  45.             s1=ii[s1p+1].split(':')[0]
  46.             s2=ii[s1p+2].split(':')[0]
  47.             if getGT(s1)==getGT(s2)!=0:
  48.                 share+=1
  49.             elif getGT(s1)!=0 and getGT(s2)==0:
  50.                 onlys1+=1
  51.             elif getGT(s1)==0 and getGT(s2)!=0:
  52.                 onlys2+=1
  53.             elif getGT(s1)==getGT(s2)==0:
  54.                 pass
  55.             else:
  56.                 diffshare+=1
  57.     print "S1 and S2 share same mutation: ",str(share)
  58.     print "only S1 mutation: ",str(onlys1)
  59.     print "only S2 mutation: ",str(onlys2)
  60.     print "S1 and S2 share diff genotype mutation: ",str(diffshare)
复制代码

  1. S1 and S2 share same mutation:  12
  2. only S1 mutation:  324
  3. only S2 mutation:  294
  4. S1 and S2 share diff genotype mutation:  4
复制代码

回复

使用道具 举报

发表回复

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则