搜索
12
返回列表 发新帖
楼主: Jimmy

从fasta文件中找出CG含量最大的DNA片段

[复制链接]

0

主题

4

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2018-5-17 14:19:09 | 显示全部楼层
我用Python写得,但不知道为什么无法读取>后的基因信息,最后一行也没显示>Rosalind_0808,而是直接显示的序列。有没有大神可以指导一下。谢谢!代码如下:
f=open("F:/GC.txt")
line=f.readline()
seq=''
header=''
for line in f:
    if line[0]=='>' and seq=='':
        header=line.strip()
    elif line[0]!='>':
        seq=seq+line.strip()
    elif line[0]=='>' and seq!='':
        seq=seq+line.strip()
        for nt in seq:
            if nt=='C' or nt=='G':
                GC=(seq.count("C")+seq.count("G"))/len(seq)
        print(header,GC)
        seq=''
header=line.strip()
GC=(seq.count("C")+seq.count("G"))/len(seq)
print(header,GC)
f.close()

运行结果如下:
0.4574468085106383
0.45918367346938777
TGGGAACCTGCGGGCAGTAGGTGGAAT 0.6091954022988506

回复 支持 反对

使用道具 举报

0

主题

5

帖子

100

积分

注册会员

Rank: 2

积分
100
QQ
发表于 2018-7-2 13:21:47 | 显示全部楼层
本帖最后由 猫寂先森 于 2018-7-2 13:28 编辑

[Python] 纯文本查看 复制代码
MAX_CG=0        # CG含量最大的DNA的CG值
MAX_title=''    # CG含量最大的DNA的名称
with open("d:/【生物信息学】/MyPrograms/ROSALIND/rosalind_GC.txt", 'r') as inputs:
    title = ''              # 用于储存循环中每次读入的DNA片段的名称
    tot = 0                 # 判断是否是第一次读入DNA名称的标志,避免计算CG值时DNA长度为0的情况出现
    CG = 0                  # 用于存储循环中每次读入的DNA片段的CG值
    line = inputs.readline()
    while  line:
        if line[0] is '>':      # 判断读入的行是否是DNA名称
            if tot == 0:        # 判断是否是第一次读入DNA名称,避免计算CG值时DNA长度为0的情况出现
                MAX_title = line[1:-1]      # 第一次读入DNA名称,直接赋值给MAX_title
                tot += 1                    # 更改标志tot的值
            else:
                CG = (DNA.count('C') + DNA.count('G')) / len(DNA) * 100     # 计算当前已读入DNA片段的CG值
                if CG > MAX_CG:
                    MAX_CG = CG             # 当当前DNA片段的CG值大于最大值时,用当前的CG值和DNA名称替换最大值
                    MAX_title = title
            title = line[1:-1]              # line读入了新的DNA名称后,替换掉当前保存的DNA的名称
            DNA = ''                        # 清空DNA序列
        elif line[-1] == '\n':              # 文件最后一行DNA序列末尾没有换行符'\n',所以需要特殊判断
            DNA += line[0:-1]
        else:
            DNA += line
        line = inputs.readline()
    else:                                   # 最后一次读入的DNA序列并未与最大值比较,所以需要再计算一次。
        CG = (DNA.count('C') + DNA.count('G')) / len(DNA) * 100
        if CG > MAX_CG:
            MAX_CG = CG
            MAX_title = title
    print(MAX_title,'\n%f'% MAX_CG)         # 输出CG值最大的DNA片段的名称和CG值


感觉python最大的障碍是读入数据太不友好了,需要判断一大堆东西。
或许是因为我才刚刚开始接触python,不太了解各种方便的包,只能自己造轮子。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

100

积分

注册会员

Rank: 2

积分
100
QQ
发表于 2018-7-2 14:08:08 | 显示全部楼层
本帖最后由 猫寂先森 于 2018-7-2 14:18 编辑
可可西里 发表于 2018-5-17 14:19
我用Python写得,但不知道为什么无法读取>后的基因信息,最后一行也没显示>Rosalind_0808,而是直接显示的 ...

不知道你用的是什么IDE,推荐一个程序的调试技巧,设置几个断点,然后用DeBug模式,用step over(中文叫步进)的方法调试,这样子可以一行一行的运行你的程序,从而可以监控到每一个变量的值的变化。
你的程序一开始line=f.readline(),这里就给line读入了一个值,是">Rosalind_****\n",但是当for循环开始的时候,又再次读了一个值给变量line,冲掉了之前的">Rosalind_****\n",所以导致输出结果时完全没有DNA片段的名称了。
然后你的第二个elif条件判断执行时又出现了 seq=seq+line.strip(),这时候的line的值可是一个">Rosalind_****\n",这将会导致你的seq以">Rosalind_****\n"结尾,seq的长度也不正确了。


下面是我用你的代码改的,结果应该是没问题了。你可以复制到自己的电脑上运行试试。
读入文件路径那里需要改成你自己的路径。
[Python] 纯文本查看 复制代码
f=open("d:/【生物信息学】/MyPrograms/ROSALIND/rosalind_GC.txt")
seq=''
header=''
for line in f:
    if line[0]=='>' and seq=='':
        header=line.strip('>\n')
    elif line[0]!='>':              # 当line[0]不为'>'时,说明此时读入line中的是DNA序列
        seq=seq+line.strip()
    else:                           # 这里只剩下line[0]=='>' and (seq != '')的情况了,其实这时候读入的line已经是下一个DNA片段的名称了
        for nt in seq:
            if nt=='C' or nt=='G':
                GC=(seq.count("C")+seq.count("G"))/len(seq)
        print(header,GC)
        header=line.strip('>\n')
        seq=''
GC=(seq.count("C")+seq.count("G"))/len(seq)
print(header,GC)
f.close()


在我电脑上的输出结果是这样子:


Rosalind_6404 0.5375
Rosalind_5959 0.5357142857142857
Rosalind_0808 0.6091954022988506

从你的代码里我也学会了strip函数的用法,感谢啦!
而且比我自己写的代码短。







回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-13 03:56 , Processed in 0.041159 second(s), 21 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.