搜索
查看: 4683|回复: 14

【菜鸟Python练习10】[ROSALIND-CONS] Consensus and Profile

[复制链接]

20

主题

68

帖子

870

积分

版主

Rank: 7Rank: 7Rank: 7

积分
870
QQ
发表于 2016-10-13 12:58:07 | 显示全部楼层 |阅读模式
本帖最后由 bioinfo.dong 于 2016-10-13 13:01 编辑

【10】Consensus and Profile
*这个题目实在不知道怎么翻译比较好
*最近忙里偷闲竟做了有10个,虽说轻车熟路,但这个过程也巩固了一些技巧。有更多的人参与就更好了。开帖的目的本是希望大家可以一起坚持把这些题目写下去,当做完最后一个题目的时候,你可以告诉别人说我已经升到Python二级菜鸟了。不管code写的好不好,都可以贴上来一起讨论,每个人写的都有闪光的地方,也都有可能给别人启示
*又到了推荐Python的时候了,有空打算专门写篇文章分析一下为什么选择Python以及如何学好Python。这次推荐的还是collections这个module,之前有个帖子里我介绍了建立有序dictionary的方法,OrderedDict,它也是collections里面的一个function。今天用到了collections里另一个function, Counter,比单纯用.count()要强大很多,一起把Python学起来吧~~






[Python] 纯文本查看 复制代码
### 10. Consensus and Profile ###

from collections import Counter

fh = open('/Users/DONG/Downloads/consesus_and_profile_output2.txt', 'wt')

rosalind = OrderedDict()
seqLength = 0
with open('/Users/DONG/Downloads/rosalind_cons.txt') as f:
    for line in f:
        line = line.rstrip()
        if line.startswith('>'):
            rosalindName = line
            rosalind[rosalindName] = ''
            continue
        rosalind[rosalindName] += line
    seqLength = len(rosalind[rosalindName])
    

A,C,G,T = [],[],[],[]
consensusSeq = ''
for i in range(seqLength):
    seq = ''
    for k in rosalind.keys():
        seq += rosalind[k][i]
    A.append(seq.count('A'))
    C.append(seq.count('C'))
    G.append(seq.count('G'))
    T.append(seq.count('T'))
    counts = Counter(seq)
    consensusSeq += counts.most_common()[0][0]


fh.write(consensusSeq+'\n')
fh.write('\n'.join(['A:\t'+'\t'.join(map(str,A)), 'C:\t'+'\t'.join(map(str,C)),
                    'G:\t'+'\t'.join(map(str,G)), 'T:\t'+'\t'.join(map(str,T))]))

fh.close()


本帖子中包含更多资源

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

x



上一篇:常见的芯片数据都可以在bioconductor里面找到注释包
下一篇:会议报告板块发帖说明
You really shouldn't spend your time reinventing the wheel
回复

使用道具 举报

1

主题

41

帖子

284

积分

中级会员

Rank: 3Rank: 3

积分
284
发表于 2016-10-14 01:23:22 | 显示全部楼层
本帖最后由 learnyoung 于 2016-10-14 01:26 编辑

试着自己用python写了一个,主要用到了numpy模块的二维数组array:
[['A' 'T' 'C' 'C' 'A' 'G' 'C' 'T']
['G' 'G' 'G' 'C' 'A' 'A' 'C' 'T']
['A' 'T' 'G' 'G' 'A' 'T' 'C' 'T']
['A' 'A' 'G' 'C' 'A' 'A' 'C' 'C']
['T' 'T' 'G' 'G' 'A' 'A' 'C' 'T']
['A' 'T' 'G' 'C' 'C' 'A' 'T' 'T']
['A' 'T' 'G' 'G' 'C' 'A' 'C' 'T']]
然后直接调出数组的每一列进行处理,版主给的Counter真的非常好用

[Python] 纯文本查看 复制代码
# coding=utf-8
import numpy as np
import os
from collections import Counter
fhand = open('dna_string1.txt')
t = []
for line in fhand:
    if line.startswith('>'):
        continue
    else:
        line = line.rstrip()
    line_list = list(line)
    t.append(line_list)
a = np.array(t)#创建一个二维数组

L1,L2,L3,L4 = [], [], [], []
comsquence=''
for i in range(a.shape[1]):
    l = [x[i] for x in a] #调出二维数组的每一列
    L1.append(l.count('A'))
    L2.append(l.count('C'))
    L3.append(l.count('T'))
    L4.append(l.count('G'))
    comsquence=comsquence+Counter(l).most_common()[0][0]
print comsquence
print 'A:',L1,'\n','T:',L2,'\n','C:',L3,'\n','G:',L4
回复 支持 2 反对 0

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-5-12 16:34:14 | 显示全部楼层
### 10. Consensus and Profile ###
按我自己的理解翻译一下:
如果只找一条链的同源,拿这条链去跟别的链一一比对,找出字符不匹配的最小的那条链就是同源链喽
但是当寻找多条链的同源的时候,这时我们需要一条所谓的“参考基准链”(就是某一位置出现概率最大的碱基组成的链,比如:给定的8条链中第一位出现概率最大的为A(出现5次),则基准链的第一位为A,第二位出现概率最大的碱基是T,则选择T,以此类推,然后将每个位置的出现概率最大的碱基连成一条链,就是参考基准链)
那个profile 是每个碱基在每条链出现的次数组成的矩阵。(可以看出A这7条链中在第一个位置出现的次数为5,第二个位置的次数为1)。应该讲明白了
我的是python3.6
[Python] 纯文本查看 复制代码
from collections import Counter
import collections

fh = open('E:/我的自学/mypyton/profile10.txt', 'wt')

rosalind=collections.OrderedDict()#这里是py3的写法,读取序列
seqLength =0
with open('E:/我的自学/mypyton/sample10.txt') as f:#提取每一个seq为line
    for line in f:
        line = line.rstrip()
        if line.startswith('>'):
            rosalindName = line
            rosalind[rosalindName] = ''
            continue
        rosalind[rosalindName] += line
    seqLength = len(rosalind[rosalindName])#每个seq的长度为8
##     
##seq是每条seq的某一位置的那个碱基连成一条seq,然后计数,排序,取最大,连成一个
A,C,G,T = [],[],[],[]#初始化
consensusSeq = ''#初始化
for i in range(seqLength):
    seq = ''
    for k in rosalind.keys():
        seq += rosalind[k]
    A.append(seq.count('A'))#计数
    C.append(seq.count('C'))
    G.append(seq.count('G'))
    T.append(seq.count('T'))
    counts = Counter(seq)#为seq计数
    consensusSeq += counts.most_common()[0][0]#从多到少返回,是个list啊,只返回第一个


    
fh.write(consensusSeq+'\n')
fh.write('\n'.join(['A:\t'+'\t'.join(map(str,A)), 'C:\t'+'\t'.join(map(str,C)),
                   'G:\t'+'\t'.join(map(str,G)), 'T:\t'+'\t'.join(map(str,T))]))#
## 
fh.close()
consensusSeq#原来输出该是这样
回复 支持 0 反对 1

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-10-13 15:43:25 | 显示全部楼层
其实我们小编团队大多数人都是用perl的~
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

20

主题

68

帖子

870

积分

版主

Rank: 7Rank: 7Rank: 7

积分
870
QQ
 楼主| 发表于 2016-10-13 23:02:37 来自手机 | 显示全部楼层
Jimmy 发表于 2016-10-13 15:43
其实我们小编团队大多数人都是用perl的~

好忧伤,哈哈,改天写篇文章拉人学Python
回复 支持 反对

使用道具 举报

20

主题

68

帖子

870

积分

版主

Rank: 7Rank: 7Rank: 7

积分
870
QQ
 楼主| 发表于 2016-10-14 02:50:08 来自手机 | 显示全部楼层
learnyoung 发表于 2016-10-14 01:23
试着自己用python写了一个,主要用到了numpy模块的二维数组array:
[['A' 'T' 'C' 'C' 'A' 'G' 'C' 'T']
[ ...

赞一个,直接用numpy的array是个非常好的idea,尤其是处理这种问题。scipy+pandas+numpy可以让Python很自如的做矩阵运算~
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-5-12 16:38:22 | 显示全部楼层
最后一句write里面的美容‘A:\t’没看懂,带我飞,求求你了
回复 支持 反对

使用道具 举报

103

主题

133

帖子

848

积分

版主

Rank: 7Rank: 7Rank: 7

积分
848
发表于 2017-8-2 14:34:20 | 显示全部楼层
本帖最后由 zckoo007 于 2017-8-2 14:37 编辑
学习最快乐 发表于 2017-5-12 16:38
最后一句write里面的美容‘A:\t’没看懂,带我飞,求求你了

\n  回车
\t   tap
.join(map(str,A))  把 A=[5, 1, 0, 0, 5, 5, 0, 0] 格式转化成 A:5 1 0 0 5 5 0 0
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-8-2 14:44:25 | 显示全部楼层
zckoo007 发表于 2017-8-2 14:34
\n  回车
\t   tap
.join(map(str,A))  把 A=[5, 1, 0, 0, 5, 5, 0, 0] 格式转化成 A:5 1 0 0 5 5 0 0 ...

回复 支持 反对

使用道具 举报

0

主题

7

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2017-11-21 23:01:27 | 显示全部楼层
附一个我的代码吧,我用的是嵌套列表,也是比较好理解的
[Python] 纯文本查看 复制代码
#*_coding: utf-8_*
def seq_list(fastqFile):
    '''读取fastq文件中的序列数据,并保存成嵌套列表'''
    seq_list = []
    for line in fastqFile.readlines():
        if not line.startswith('>'):
            seq = list(line.strip())
            seq_list.append(seq)
    return seq_list

def statistic_base(seq_list):
    for base_type in 'AGCT':
        base_total = []
        for sit in range(len(seq_list[0])):
            column = [x[sit] for x in seq_list]
            number = column.count(base_type)
            base_total.append(number)
        print '%s:%s'%(base_type, base_total)


fastqFile = open('E:\\bioinfo\data\consensusu_and_profile.fastq')
seq_list = seq_list(fastqFile)
statistic_base(seq_list)
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-26 19:04 , Processed in 0.047093 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.