搜索
查看: 2202|回复: 0

生信编程实战代码注释【第二题-hg19基因组序列的一些探究】

[复制链接]

10

主题

11

帖子

400

积分

中级会员

Rank: 3Rank: 3

积分
400
发表于 2017-7-2 00:14:07 | 显示全部楼层 |阅读模式
本帖最后由 默存 于 2017-7-2 00:15 编辑

[Python] 纯文本查看 复制代码
'''
以下代码来自【生信编程直播第二题-hg19基因组序列的一些探究-python-李治鑫】如发现错误,请在后面指出,以免祸害后来者
注释者:默存
20170701
'''

# coding: utf-8
import os 
os.chdir("E:\jupyter\文件\第二题")
from collections import OrderedDict
chr_dict = OrderedDict()
temp_chr = ""

with open("test.fasta","r") as f:
    for line in f:
        line = line.rstrip() 
        if line.startswith(">"):#如果行头是>,此行为非序列行
            tem_chr = line#把>行存储在tem_chr中
            chr_dict[tem_chr] = ""#把tem_chr导入到字典里
        else:
            chr_dict[tem_chr] += line#迭加各序列行,得完整序列
#print(chr_dict),可以print下,看看chr_dict是什么样子的

for seqName,seq in chr_dict.items():# 使用字典items()函数,取出键值对
    seqLen = len(seq)  #len()函数求序列总长
    A_count = seq.count("A") + seq.count("a") #count()函数算碱基个数,注意每种碱基可能有大小写
    T_count = seq.count("T") + seq.count("t")
    C_count = seq.count("C") + seq.count("c")
    G_count = seq.count("G") + seq.count("g")
    N_count = seq.count("N") + seq.count("n")
    GC = (C_count + G_count)/(seqLen - N_count)#计算GC比例,一般忽略掉未知碱基N
    
    print("A_count:","%2.f"%A_count,"T_count:","%2.f"%T_count,"G_count:","%2.f"%G_count,"C_count:","%2.f"%C_count\
    ,"N_count:","%2.f"%N_count,"%GC:","%5.2f"%(GC),"seqLen:","%2.f"%(seqLen))#格式化输出,让输出整洁点,"%5.2f"5表示字宽,2表示保留2位小数
    
    
'''
以下是使用python中的pysam模块完成,需要先安装pysam模块,代码来自【生信编程直播第二题-hg19基因组序列的一些探究-python-李治鑫】
注释者:默存
20170701
'''
# -*- coding: UTF-8 -*-  
import pysam
import os
#os.chdir("./")
f = pysam.FastaFile("test.fasta")#pysam提供的方法打开文件 
#dir(hg19) #list methods
#list(zip(hg19.references,hg19.lengths))#zip()迭代多个列表或者元祖
for seqName in f.references:
        seq = f.fetch(seqName)#fetch取出对应名字的序列
        #字符串len()函数求序列总长
        seqLen = len(seq)  
        A_count = seq.count("A") + seq.count("a") #count()函数算碱基个数,注意每种碱基可能有大小写
        T_count = seq.count("T") + seq.count("t")
        C_count = seq.count("C") + seq.count("c")
        G_count = seq.count("G") + seq.count("g")
        N_count = seq.count("N") + seq.count("n")
        GC = (C_count + G_count)/(seqLen - N_count)
        
        print("A_count:","%2.f"%A_count,"T_count:","%2.f"%T_count,"G_count:","%2.f"%G_count,"C_count:","%2.f"%C_count\
    ,"N_count:","%2.f"%N_count,"%GC:","%5.2f"%(GC),"seqLen:","%2.f"%(seqLen))
    
'''
这一节补充学习的:
fasta数据格式
又熟悉了一次字典中的dict.keys(),dict.values(),dict.items()
字符串的count()函数,len()函数
格式化输出
pysam模块下周真的要转行去做生信分析工作了,好期待,好怕怕!
'''    






上一篇:【免疫组测序-2】免疫组测序背景介绍(免疫组库)
下一篇:落入窠(ke)臼(jiu):GATK best practice每个步骤都是必须的吗?
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-20 13:28 , Processed in 0.029245 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.