![e03f8aafa643b3507b6e2730d72d306d.png](https://i-blog.csdnimg.cn/blog_migrate/96e874814f0a3ce629b1885877b24b3d.jpeg)
前言:Hello,各位老铁,经过上一个帖子《生物信息学中必须掌握的拍照技巧》一文,相信大家对我应该有所了解,本人本着蹭流量,蹭热度,不负责的态度,又来无耻地骗取大家的流量了~(当你点开这篇文章的时候,本文阅览量又+1)。
事先声明,本文属于小白向的文章,各位编程已经熟练的大佬们可以点返回键退出本文了~
正文:
当我们日常在跑流程的时候,我们经常会遇到软件输入文件格式与手头上的文件格式不符,面对着少则几千,多着几万(单个文件)的序列信息,顿觉无泪,又不能愉快的玩耍了。
这个时候,我们就必须请出我们有力的生产工具----Python字典,对那几千几万条信息进行整理。
对于不明白Python字典的童鞋,我们可以在这里复习一下,字典是个什么东西?
字典(Dictionary)在Python中是一种可变的容器模型,它是通过一组键(key)值(value)对组成,这种结构类型通常也被称为映射,或者叫关联数组,也有叫哈希表的。每个key-value之间用“:”隔开,每组用“,”分割,整个字典用“{}”括起来,用看得懂的代码表示就像下面这样:
dict_1 = { ‘隔壁老侯’: ‘帅哥’, ’契卡’ :’催稿狂’, ’小孟孟’ :’天然萌’ }#直接创建
又或者是这样:
items=[(’隔壁老侯’ , ‘帅哥’), (’契卡’, ‘催稿狂’),( ’小孟孟’, ‘天然萌’)]
dict_2 = dict(items)#利用dict创建
再或者是这样:
dict3 = dict(’隔壁老侯’ = ‘帅哥’, ’契卡’ = ‘催稿狂’, ’小孟孟’ = ‘天然萌’)#另一种dict创建方法
怎么去理解这堆东西呢?
其实很简单,映射,就像我们小学语文的一对一连线题,一个key对应一个value,就像提到隔壁老侯,我们会不由自主的想起他是一个帅哥,而想到契卡,我们会想到他在群里每天催稿,而提起孟大哥,我们总是情不自禁的想到他卖萌的瞬间,像北大浩巍,在线发嗲之类的经典片段。这么说应该会比较理解了吧~
接下来我们来看看字典的几大常规操作,
福斯特汪~(来自单身狗的咆哮~),怎么输出字典里的值,很简单,把key输出来就是了,怎么理解,再拿上面的字典来解释,这个时候,我们添加入一个情景,大家跟着我一起来进入语文大课堂:
语文课上,老师提问:“隔壁老侯是怎么样的人啊~?”
学生回答:“帅哥~~~~”
就是这么简单明了,具体是怎么样的骚操作呢,如下所示:
dict_1 = { ‘隔壁老侯’: ‘帅哥’, ’契卡’ :’催稿狂’, ’小孟孟’ :’天然萌’ }
Print(‘隔壁老侯是一个’,dict[‘隔壁老侯’])#就相当于在字典dict1里面提取“隔壁老侯”的索引(此时学生的脑子里飘过一个帅哥的形象)
然后在解释器里面就出现下面的字眼:
隔壁老侯是一个帅哥
接着,我们来谈一下字典的增减元素,这个就像班里转来一名同学或者转走一名同学,究竟怎么操作,具体如下:
dict_1 = { ‘隔壁老侯’: ‘帅哥’, ’契卡’ :’催稿狂’, ’小孟孟’ :’天然萌’ }
dict1[‘二货潜’] = ‘搬运小王子’
此时如果我们print一下dict1,那么就会得到如下所示的效果:
dict_1 = { ‘隔壁老侯’: ‘帅哥’, ’契卡’ :’催稿狂’, ’小孟孟’ :’天然萌’, ‘二货潜’=’搬运小王子’ }
有时候我们发现自己添加的东西错了,然后要修改一下,就像我们发现孟大哥的才华已经把他天然萌的属性给遮盖了,这个时候,我们就需要对字典进行修改,
dict1[‘小孟孟’] = ‘大神’
这个时候,我们再print一下字典,就会发现之前要修改的东西已经修改了,
dict_1 = { ‘隔壁老侯’: ‘帅哥’, ’契卡’ :’催稿狂’, ’小孟孟’ :’大神’, ‘二货潜’=’搬运小王子’ }
然后我们又发现,隔壁老侯除了帅就一无是处了,他不能跟这些人在同一个字典里面,我们要把这种金玉其外败絮其中的禽兽剔除掉,这个时候我们就需要删除字典里面的元素,具体操作如下:
del dict1[‘隔壁老侯’]
然后我们再print一下,就会发现这个衣冠禽兽不见了
dict_1 = { ’契卡’ :’催稿狂’, ’小孟孟’ :’大神’, ‘二货潜’=’搬运小王子’ }
然后,这几个人在毕业以后全部成功跳坑成功,不再是一头可怜的生物狗,亦或者生信猿,这个时候,我们的字典里面就没什么东西了,操作如下
dict1.clear()
如果再print这个字典,我们会发现we have anything(你个穷*)~,
那么一个空的字典还有什么用,删了删了,留着占内存,
del dict1
如果我们再print,我们会发现它已经消失了~
下面还有一些常用的用法,你们就看看吧,偶尔会用得上:
items()#获取字典键值对数据,以列表形式返回;
keys()#获取字典键的数据,以列表形式返回;
values()#获取字典值的数据,以列表形式返回;
复习完基础知识,那么我们就进入今天的主体吧,如何利用字典进行文件处理
事情是这样的,最近我在做物种亲缘鉴定的时候用到了blat(不是blast,究竟是什么可以自行百度一下),得到的结果图如下。
![a41fe38fd2db05a0a955fc60055c73ae.png](https://i-blog.csdnimg.cn/blog_migrate/e3ee0c82489c4d818330777ff158c193.jpeg)
而我要从这张表里面提取信息,就是我圈出来的地方,
首先就是q-name,是我blat比对时作为sample的序列的unigene的序列名,由于一条序列可能比对上多条reference上面的序列,所以只能根据其他提供的信息进行过滤,得到最佳的序列,
其次就是t-name,是blat时作为reference的unigene的序列名,
接着就是match跟mis-match,数值都代表着碱基数,
另一个重要的参数blockcount则是根据它们mapping的时候出现的gap数量得到的一个值,blockcount为1表示没有gap的情况,随着gap增加blockcount的数值会增大。
基本信息就如下,那么接下来就是我们的伪代码阶段,理清一下我们需要怎么写代码:
首先,我们知道q-name可能出现相同重复的序列,而我们要得到唯一的,那么最好是建立一个一对一的映射,然后在出现相同的序列名时进行判断即可,那么此时,我们是的初步代码可以这么写:
import sys#导进sys库
x = sys.argv[1]#选择传入第二个参数(0是第一个,故而1是第二个),即我们在运行代码的时候在命令行里面这样写 python 脚本.py(第一个参数) 输入文件(第二个参数)> 结果
file = open(x, 'r')#以可读的形式读取传入参数所代表的文件
lines = file.read()#以read()的形式读入
dic = {}#创建一个空字典
keys = []#创建一个list存储key,方便后面判断的时候使用
for line in lines.split('n'):#将读进来的东西以换行符为分割,此时效果相当于readlines,只不过一次性读入内存
if q-name第一次出现:
q-name作为key,然后其他的作为value写入
出现的q-name也写入keys里面,方便以后判断用
if q-name再次出现:
if 是更加合适的序列:
更新value的内容(不是更合适的就不更新,就不用写了)
print 结果
接着出现在我们面前的一个问题就是,这个更加合适的序列要怎么判断?
首先,blockcount为1代表着没有出现gap,一条序列如果能完美契合上去,证明是完全一致的;如果gap多的话,那mis-match肯定也随着增多,而整条序列的长度如果特别长,那出现几个mis-match显然对blockcount的影响会比较小,而长度较短,则影响大一些,match跟mis-match的优先级显然低一些,所以优先判断blockcount,然后再判断match跟mis-match。
接着,blockcount要怎么判断,由表我们可以知道(虽然我在前面也说了,但不妨仔细观察一下),然后gap跟mis-match的数量增多,那blockcount的数值就会增大,而1代表着完美的匹配,那么我们就知道,优先级肯定是由小到大,如果出现了1,那其他数值就不用看了。
然后的话就是match跟mis-match,它们代表着碱基数量的多少,当blockcount相同时,那还就需要match跟mis-match来辅助判断,很明显是碱基数越多的两条序列会更接近一些,所以match跟mis-match的总数大的优先级高一些。
好了,基本的判断就算完了,接着,我们跟着这个思路,把伪代码写出来吧
import sys
x = sys.argv[1]
file = open(x, 'r')
lines = file.read()
dic = {}
keys = []
for line in lines.split('n'):
if q-name第一次出现:
q-name作为key,然后其他的作为value写入
出现的q-name也写入keys里面,方便以后判断用
if q-name再次出现:
if blockcount比之前的小:
更新value的内容
elif blockcount跟之前的一样:
if match跟mis-match的总和比较大:
更新value的内容
print 结果
基本这样,我们的伪代码就算完成了,接着是如何把他转换为可以执行的python代码,我们还是来看一下这张图:
![9a74b326e8fd96a5e2e1190f8ac7eaaa.png](https://i-blog.csdnimg.cn/blog_migrate/59b0096ba1fa50440d68de075641017b.jpeg)
首先,我们数一下我们需要提取的东西的位置,可以看到match在第一位,在代码里面我们在怎么实现,
记得我们前面只是将整个文件读进来,然后分割成一行一行的,那么我们现在要做的就是把一行里面分割成数个模块,整一行其实是以tab键进行分割的(不信的可以用UE打开看看),所以我们对整一行进行分割就要以tab键进行分割,即:
line.split('t')
这样我们就把一行以tab键分割成了n段,至于有多少呢,数数tab键吧,一行有21个,而我们所需要的分别是0,1,9,13,17,其中的9作为字典的key,而其余都写入value,看到这里,可能有些童鞋懵逼了,怎么写???
回想一下,我们刚刚的基础知识
dict[key]=value
而我们现在其实就是
dict[q-name]=t-name+match+mis-match+blockcount,
而在代码里面怎么实现,接着看就知道了,
我们上面将一行分割为n列,为了方便我们后面写,一般的话,我们可以在分割的时候,令某个变量= line.split('t'),这里的话,我习惯用strings,即写成这样:
strings=line.split(‘t’)
而对应的q-name就可以写成strings[9]#即strings这里面,以tab键分割的第十个元素,而其他的也可以依次写出来,即t-name为strings[13],match为strings[0]等等。
接着,我们把代码块完善一下吧,
import sys
x = sys.argv[1]
file = open(x, 'r')
lines = file.read()
dic = {}
keys = []
for line in lines.split('n'):
strings = line.split(‘t’)
if strings[9] not in keys:#strings[9]即q-name第一次出现
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
keys.append(strings[9])
if strings[9] in keys:
if strings[17] < int(dic[strings[9]]).split(‘>’)[3]:#将dic[strings[9]]对应的value以>>分割,然后将第四个元素(注意观察前面的字典的value里面blockcount位于第四)转化为整数型跟新出现的blockcount比较
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
elif strings[17] = int(dic[strings[9]]).split(‘>’)[3]:
if int(strings[0]) + int[strings[1]]> int(dic[strings[9]]).split[1] + int(strings[9]).split[2]:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
print 结果
好了,基本整个代码的核心框架就算完成,接着来抠一下细节,首先的话,先判断结构就写得不够简洁,可以把第一次出现的写成else,这样更加整洁舒服一些,
import sys
x = sys.argv[1]
file = open(x, 'r')
lines = file.read()
dic = {}
keys = []
for line in lines.split('n'):
strings = line.split(‘t’)
if strings[9] in keys:
if strings[17] < int(dic[strings[9]]).split(‘>’)[3]:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
elif strings[17] = int(dic[strings[9]]).split(‘>’)[3]:
if int(strings[0]) + int[strings[1]]> int(dic[strings[9]]).split[1] + int(strings[9]).split[2]:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
else:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
keys.append(strings[9])
接下来,我们就可以考虑一些打印出来的效果会怎么样,如果直接打印的话,那就会出现{‘key’:’value’,’key’:’value’}这样,那么我们就需要观察,然后把这些杂七杂八的符号给替换掉,最后的代码放给大家,大家可以在理解的基础上,找一些诸如blast的结果表格的练习一下,最后的代码如下(请需要练习的各位一个字一个字的打出来,代码要自己熟练了才有用,不然全是别人的):
import sys
x = sys.argv[1]
file = open(x, 'r')
lines = file.read()
dic = {}
keys = []
for line in lines.split('n'):
if "TRINITY" in line:#如果不判断的话,会提示index out of range,因为文件的每一行分割并不是一样的
strings = line.split(‘t’)
if strings[9] in keys:
if strings[17] < int(dic[strings[9]]).split(‘>’)[3]:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
elif strings[17] = int(dic[strings[9]]).split(‘>’)[3]:
if int(strings[0]) + int[strings[1]]> int(dic[strings[9]]).split[1] + int(strings[9]).split[2]:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
else:
dic[strings[9]] = strings[13] + ‘>’ +strings[0] + ‘>’ + strings[1] + ‘>’ + strings[17]
keys.append(strings[9])
rlt = str(dic)
rlt = rlt.replace("{'", '')
rlt = rlt.replace("': '", "t")
rlt = rlt.replace(">", "t")#如果你在前面字典的分隔符就用t,等你最后出来的结果t不会变成tab键,而还是t的分隔符,故而在字典里面推荐使用其他分隔符,方便后面转换
rlt = rlt.replace("', '", "n")
rlt = rlt.replace("'}","n").rstrip()
print("Qname" + "t" + "Tname" + "t" + "match" + "t" + "mis-match" +"t" + "blockSizes")
print(rlt)
![1cd6c6d8b1a9b724f8b9e6d9f004df75.png](https://i-blog.csdnimg.cn/blog_migrate/73e98382d5a95e370d09e29a7c6e65c1.jpeg)
如果替换完成,最后的结果会比较整洁,像这样:
![f6f1ff1d8c8b2cba76c880926ad5b593.png](https://i-blog.csdnimg.cn/blog_migrate/6c44255759f9bc7a7d3f5e0d47e33933.jpeg)
后记:其实使用字典对这些多条件的文件进行过滤的效率并不好,特别是if语句的使用会更加慢,如果对pandas比较熟悉的话,用pandas写会比较有效率,我也要努力学习pandas,争取早日写一篇关于pandas在我们生信日常使用中的应用~就这样啦,我们下回再见~
练习文件见此:链接:https://pan.baidu.com/s/1Vm3MK_u-89F7jFzuV9JE0g
提取码:ow32