生信编程实战第1题(python)

题目来自生信技能树
统计人类外显子长度
坐标的文件可如下下载：

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt

打开文件如下

image.png

于是写了如下脚本

import sys
import re
args=sys.argv
filename=args[1]
exon_length=0
aDict={}
with open (filename) as fh :
      for line in fh:
          if line.startswith("#"):
              continue
          lineL=line.strip().split("\t")
          exon_position=lineL[-2] #取出倒数第二列，坐标列
          if exon_position=="-":  #有的基因没有外显子的坐标，用-代替的，所以这行就要除掉，不然会报错
              continue
          exon_position=re.sub("\[|\]","",exon_position) #把坐标列的[]去除，注意正则表达式的用法
          exonL=exon_position.split(",")  
          for exon in exonL:
              exonS=lineL[0]+":"+exon    #有点基因会有相同坐标的外显子，所以要去除这一部分，注意要加上染色体的编号，染色体不同而坐标一样就没事
              if exonS not in aDict:    #如果坐标没有在字典，即第一次出现，就将其放入字典，并继续操作。
                   aDict[exonS]=1
                   exon_pL=exon.split("-")
                   exon_start=int(exon_pL[0].strip())
                   exon_end=int(exon_pL[1].strip())
                   exon_length+=exon_end-exon_start
print(exon_length)

然后运行

ubuntu@VM-0-4-ubuntu:~/data/practice$ python3 exon.py CCDS.current.txt 
36443621

正则表达式再学习熟悉熟悉