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

bioinfo.dong

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






[mw_shl_code=python,true]### 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]
    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()[/mw_shl_code]

learnyoung

试着自己用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真的非常好用

[mw_shl_code=python,true]# 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 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
[/mw_shl_code]

学习最快乐

### 10. Consensus and Profile ###
按我自己的理解翻译一下：
如果只找一条链的同源，拿这条链去跟别的链一一比对，找出字符不匹配的最小的那条链就是同源链喽
但是当寻找多条链的同源的时候，这时我们需要一条所谓的“参考基准链”（就是某一位置出现概率最大的碱基组成的链，比如：给定的8条链中第一位出现概率最大的为A（出现5次），则基准链的第一位为A，第二位出现概率最大的碱基是T,则选择T，以此类推，然后将每个位置的出现概率最大的碱基连成一条链，就是参考基准链）
那个profile 是每个碱基在每条链出现的次数组成的矩阵。（可以看出A这7条链中在第一个位置出现的次数为5，第二个位置的次数为1）。应该讲明白了
我的是python3.6
[mw_shl_code=python,true]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#原来输出该是这样[/mw_shl_code]

Jimmy

其实我们小编团队大多数人都是用perl的~

bioinfo.dong

Jimmy 发表于 2016-10-13 15:43
其实我们小编团队大多数人都是用perl的~

好忧伤，哈哈，改天写篇文章拉人学Python

bioinfo.dong

learnyoung 发表于 2016-10-14 01:23
试着自己用python写了一个，主要用到了numpy模块的二维数组array：
[['A' 'T' 'C' 'C' 'A' 'G' 'C' 'T']
[ ...

赞一个，直接用numpy的array是个非常好的idea，尤其是处理这种问题。scipy+pandas+numpy可以让Python很自如的做矩阵运算～

学习最快乐

最后一句write里面的美容‘A:\t’没看懂，带我飞，求求你了

zckoo007

学习最快乐 发表于 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

学习最快乐

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 ...

biostudy

附一个我的代码吧，我用的是嵌套列表，也是比较好理解的
[mw_shl_code=python,true]#*_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)[/mw_shl_code]