搜索
楼主: Jimmy

生信编程直播第0题-生信编程很简单!

[复制链接]

0

主题

16

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-15 12:46:54 | 显示全部楼层
本帖最后由 azazelcc 于 2017-5-20 14:37 编辑

有一种感觉到期末,狂补作业的感觉。
因为被参数处理包的面对对象搞蒙了,用了一个月时间学习了一下python的面对对象,现在看到“.”神清气爽了。

以下是python3代码

大部分都是按照dongye老师的方法,其中长度分布用matplotlib包做了直方图。
Update1:发现bug,因为fasta的demo文件中有N缺失碱基,原程序并未对此做出过滤,会报错。(待改正)Update2:其实处理缺失碱基还是个挺有讲究的东西,我这里就为了偷懒,直接remove了N碱基,希望这点可以跟大家讨论。

[Python] 纯文本查看 复制代码
#!/usr/bin/env python
# encoding=utf-8

"""
My first bioinformatic tools for FASTA file

Richard 2017
"""

from optparse import OptionParser

import sys

args = sys.argv


def readfasta(filename):
"""
:param filename: ;要读取的FASTA文件名
:return: 返回基因名:序列的字典文件
"""

fa = open(filename, 'r')
res = {}
ID = ''
for line in fa:
if line.startswith('<'):
ID = line.strip('\n').strip('N')
res[ID] = ''
else:
res[ID] += line.strip('\n')

return res


def cut35(filename, cut_5, cut_3, fileout):
"""
用于从5/'和3/'端分别去掉几个碱基
:param filename: 读入的文件
:param cut_5: 5/'端要去掉的碱基数量
:param cut_3: 3/'端要去掉的碱基数量
:return: 返回文件fileout
"""

of = open(fileout, 'w')

data = readfasta(filename)

for key, value in data.items():
of.write(key + '\n')
out = value[int(cut_5):][:-int(cut_3)]
of.write(out + '\n')


def distribution(filename):
"""
用直方图表示序列长度的分布
用到numpy和matplotlib包
:param filename: 
:return: 
"""

import numpy as np
import pylab as pl

data = readfasta(filename)
res = []

for value in data.values():
res.append(len(value))

dataarray = np.array(res)

pl.hist(dataarray)
pl.xlabel('The length of sequences')
pl.show()


def countGC(filename, fileout):
"""
统计整个FASTA文件中的GC含量和占比
:param filename: 
:param fileout: 
:return: 
"""
fo = open(fileout, 'w')
data = readfasta(filename)
NUM = 0
GContent = 0
for key, value in data.items():
value = value.upper()
NUM += len(value)
GContent += (value.count('G') + value.count('C'))
output = 'The total length is ' + str(NUM) + '\n' + 'The total GC content is ' + str(
GContent) + '\n' + 'The GC ration is ' + str(float(GContent / NUM))
print(output)
fo.write(output)


def complement(filename, fileout):
"""
取互补序列,即DNA to RNA
:param filename: 
:param fileout: 
:return: 
"""

fo = open(fileout, ' w')
data = readfasta(filename)
seq = ''

compdic = {'A': 'A', 'T': 'U', 'G': 'G', 'C': 'C'}

for key, value in data.items():
fo.write(key + '\n')
for i in value:
seq += compdic[i]
fo.write(seq + '\n')


def reverse(filename, fileout):
"""
取反向序列
:param filename: 
:param fileout: 
:return: 
"""
fo = open(fileout, 'w')
data = readfasta(filename)
seq = []

for key, value in data.items():
fo.write('%s(reverse)\n' % key)
seq = list(value)
seq = seq.reverse()
seq = ''.join(seq)
fo.write('%s\n' % seq)


def main(args):
parser = OptionParser()
parser.add_option("-f", "--fasta", dest="filename", help="fasta filename",
metavar="FILE")
parser.add_option("-o", "--txt", dest="fileout", help="output filename",
metavar="FILE")
parser.add_option("--cut", "--cut-fasta", dest="cut", help="cut fasta file",
action="store_true", default=False)
parser.add_option("-5", "--cut-5", dest="cut_5", type=int, help="cut fasta N(N>0) from 5'",
metavar="INT", default=0)
parser.add_option("-3", "--cut-3", dest="cut_3", type=int, help="cut fasta N(N>0) from 3'",
metavar="INT", default=0)
parser.add_option("--len-dis", dest="len_dis", help="sequence length distribution",
action="store_true", default=False)
parser.add_option("--gc", "--count-gc", dest="gc", help="pring GC%", action="store_true",
default=False)
parser.add_option("--comp", "--complement-seq", dest="comp", help="get complement sequences", action="store_true",
default=False)
parser.add_option("--rev", "--reverse-seq", dest="rev", help="get reverse sequence", action="store_true",
default=False)

(options, args) = parser.parse_args()

if not options.filename:
parser.print_help()
filename = options.filename
fileout = options.fileout

if options.cut:
cut_5 = options.cut_5
cut_3 = options.cut_3
cut35(filename=filename, cut_5=cut_5, cut_3=cut_3, fileout=fileout)

if options.len_dis:
distribution(filename=filename)

if options.gc:
countGC(filename=filename, fileout=fileout)

if options.comp:
complement(filename=filename, fileout=fileout)

if options.rev:
complement(filename=filename, fileout=fileout)


if __name__ == '__main__':
main(args)






回复 支持 反对

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2017-5-15 20:15:59 | 显示全部楼层
azazelcc 发表于 2017-5-15 12:46
有一种感觉到期末,狂补作业的感觉。
因为被参数处理包的面对对象搞蒙了,用了一个月时间学习了一下python ...

这么长的代码呀
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

16

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-15 23:29:12 | 显示全部楼层
Jimmy 发表于 2017-5-15 20:15
这么长的代码呀

写了一些注释,然后有个坏习惯就是爱换行。还有最后一题,作业就全补上了!
回复 支持 反对

使用道具 举报

16

主题

61

帖子

515

积分

高级会员

Rank: 4

积分
515
发表于 2017-6-10 06:45:48 | 显示全部楼层
anlan 发表于 2017-3-27 21:20
034-perl-shell

输入文件:longreads.fq

第40行应该是max
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-6-11 13:35:01 | 显示全部楼层

嗯,打错了,好细心
回复 支持 反对

使用道具 举报

16

主题

61

帖子

515

积分

高级会员

Rank: 4

积分
515
发表于 2017-6-11 17:42:17 | 显示全部楼层
anlan 发表于 2017-6-11 13:35
嗯,打错了,好细心

试验了一下,并改造了一下您的代码,看到的,希望多多与您交流!也希望您能多多帮忙注释代码,促进学习。谢谢!
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-6-11 18:43:55 | 显示全部楼层
泥土 发表于 2017-6-11 17:42
试验了一下,并改造了一下您的代码,看到的,希望多多与您交流!也希望您能多多帮忙注释代码,促进学习。 ...

顶,加油
回复 支持 反对

使用道具 举报

11

主题

52

帖子

279

积分

中级会员

Rank: 3Rank: 3

积分
279
发表于 2017-6-15 11:11:56 | 显示全部楼层
xxnku 发表于 2017-5-13 16:59
老师你好,第43行如果写成print line_2[cut_5:-cut_3] 有没有什么不妥

没什么不同的
回复 支持 反对

使用道具 举报

0

主题

3

帖子

33

积分

新手上路

Rank: 1

积分
33
发表于 2017-7-27 16:43:48 | 显示全部楼层
来学习了,大家的代码写的很好啊
回复 支持 反对

使用道具 举报

0

主题

2

帖子

39

积分

新手上路

Rank: 1

积分
39
发表于 2017-10-16 16:35:56 | 显示全部楼层
叹为观止了,感觉学起来挺不容易的,入门容易,入了门就全是大深坑了,一定要加油!!!
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 15:38 , Processed in 0.039835 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.